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Aaron Christopher Boley 
THE THREE-DIMENSIONAL BEHAVIOR OF 
SPIRAL SHOCKS IN PROTOPLANETARY DISKS 

ABSTRACT: In this dissertation, I describe theoretical and numerical studies 
that address the three-dimensional behavior of spiral shocks in protoplanetary disks 
and the controversial topic of gas giant formation by disk instability. For this work, I 
discuss characteristics of gravitational instabilities (GIs) in bursting and asymptotic 
phase disks; outline a theory for the three-dimensional structure of spiral shocks, 
called shock bores, for isothermal and adiabatic gases; consider convection as a source 
of cooling for protoplanetary disks; investigate the effects of opacity on disk cooling; 
use multiple analyses to test for disk stability against fragmentation; test the sensi- 
tivity of CI behavior to radiation boundary conditions; measure shock strengths and 
frequencies in Gl-bursting disks; evaluate temperature fluctuations in unstable disks; 
and investigate whether spiral shocks can form chondrules when GIs activate. The 
numerical methods developed for these studies are discussed, including a radiation 
transport routine that explicitly couples the low and high optical depth regimes and 
a routine that models ortho and parahydrogen. Finally, I explore the hypothesis that 
chondrule formation and the FU Ori phenomenon are driven by GI activation in dead 
zones. 
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Chapter 1 



INTRODUCTION 

The formation of planetary systems in disks was anticipated as early as 1755 by 
Kant. Since that time, analytic and numerical work has demonstrated that disks 
should form during protostellar core collapse (e.g., Ulrich 1976; Cassen & Moosman 
1981; Durisen et al. 1989; Yorke et al. 1993; Pickett et al. 1997; Vorobyov & Basu 
2006; Krumholz et al. 2007), and direct imaging and spectral energy distribution 
(SED) fitting have verified that disks do surround young stars (e.g., Beckwith et al. 
1990; Strom et al. 1993; Padgett et al. 1999; Calvet et al. 2005; Andrews & WiUiams 
2005; D'Alessio et al. 2006; Eisner & Carpenter 2006). It is in these disks where grain 
growth and annealing occur (Natta et al. 2007), where planetesimals form from the 
processed dust, and where planets are built. Therefore, these protoplanetary disks 
do not simply act as a mass reservoir for accretion onto the star, but also serve 
as an astrophysical factory for chemical and dynamical processes leading to planet 
formation. For this discussion, I will focus on T Tauri-mass stars and their disks, 
although Herbig Ae/Be stellar disks are likely to be susceptible to similar phenomena. 

1 
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T Tauri stars are pre-main-sequence (PMS) stars that lie between the birth- 
hne, where PMS stars first begin quasi-static gravitational contraction (Stabler 1983, 
1988), and the zero-age main sequence (ZAMS), and have masses < 2.5Mq (Lawson et 
al. 1996). These objects are typically divided into two groups: classical T Tauri stars 
(CTTS) and weak-line (weak-emission) T Tauri stars (WTTS). CTTSs have excess 
infrared emission, excess blue continuum emission, and strong line emission. These 
data are interpreted to indicate the presence of a disk and active accretion onto the 
star. WTTSs have some infrared excess, but much weaker hne and blue continuum 
emission; the gaseous inner disk surrounding the star is believed to be mostly dissi- 
pated. The distinction between a CTTS and a WTTS is quantified by the equivalent 
width (EW) of the Ha line, where WTTSs have an EW(Hq;) < 10 A (see Andrews 
& Carpenter 2005 for a summary). The higher mass analogs of T Tauri stars are the 
Herbig Ae/Be stars (Herbig 1960). These lie between the intermediate-mass birth- 
hne (Palla & Stabler 1990) and the ZAMS, exhibit infrared excess, and have strong 
emission lines (The et al. 1994; van den Ancker et al. 1997). 

T Tauri stars fit into a more general classification of PMS objects, the Young 
Stellar Objects (YSOs). The term YSO should be understood to include the forming 
star, disk, and envelope, when present. Typically, YSOs are divided into three general 
classifications (Lada & Wilking 1984; Adams et al. 1987; Kenyon & Hartmann 1987; 
Lada 1987; Greene et al. 1994; Andrews & Williams 2005): Class I objects are embed- 
ded protostars with envelope material that is accreting onto a surrounding disk. For 



3 

these sources, energy dissipation by disk accretion is an important component of the 
total disk emission. In Class II objects, most of the disk emission is reprocessed star 
light, and envelope accretion has essentially ceased. Finally, Class III objects have 
very little disk emission. In addition to Classes I through III, Andre et al. (1993) 
proposed the categorization Class 0, with VLA 1623 as the prototype. They suggest 
that Class sources are the precursors to the Class I stage and that most of the 
circumstellar material is distributed in an envelope. In contrast, Jayawardhana et 
al. (2001) argue that Class and Class I YSOs may be at a similar stage in evolution, 
with Class YSOs forming in a denser environment than Class I YSOs. 

Classes I, II, and III can be specified quantitatively by calculating the spectral 
index, a — dlni/F^/dlni/, from the object's SED. Following Greene et al. (1994), 
Class I objects have a negative a as measured between 2.2 and 10 /xm. Class II 
objects have a flat or rising spectrum with < a < 1.6, and Classes III have an 
a > 1.6. CTTSs span Classes I and II, and WTTSs are Class III objects. The 
classification is interpreted as an evolutionary sequence, but absolute time intervals 
are not necessarily associated with a given phase. Recent observations by Eisner & 
Carpenter (2006) seem to indicate that a Class I object evolves to a Class II object 
in about 1 Myr. 

It is generally accepted that during the Class I stage of T Tauri star evolution, 
disks are massive enough for gravitational instabilities (GIs) to develop (see below). 
For example, recent observations of FU Ori, which is transitioning from the Class 
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I to II stage, indicate that the inner 1 AU may be gravitationally unstable (Zhu et 
al. 2007). Moreover, detailed measurements of L1551 IRS 5, the prototypical Class 
I object (Adams et al. 1987), indicate that the circumstellar disk mass around the 
northern source is comparable with the star's mass and that it too is gravitationally 
unstable (Osorio et al. 2003). Three-dimensional hydrodynamics simulations by Cai 
et al. (2007) of L1551 based on the disk parameters reported by Osorio et al. confirm 
that GIs develop in the disk. 

A gravitational instability is a dynamic instability driven by self-gravity. Toomre 
(1964) demonstrated that an infinitely thin gaseous disk is unstable to self-gravitating 
ring instabilities when the Toomre Q parameter. 



approaches unity. The stabilizing quantities are the sound speed Cg and the epicychc 
frequency of the gas k, and the destabilizing quantity is the surface density E. When 
Q < 1.7, a disk with finite thickness is susceptible to nonaxisymmetric instabilities 
(see Durisen et al. 2007a for a review), and spiral waves develop. For a Gl-active disk, 
the spiral waves driven by self-gravity may be the dominate way angular momentum 
is transfered outward and mass transfered inward (Lynden-Bell & Kalnajs 1972; Boss 
1984b; Larson 1984; Durisen et al. 1986; Laughlin & Bodenheimer 1994). However, 
the role that GIs play in T Tauri disks of ages > 1 Myr is uncertain, especially their 
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role in planet formation. Are disks that are older than 1 Myr cold and/or massive 
enough for GIs to activate? 

Even in relatively low-mass disks, GIs may activate for large r. A back-of-the- 
envelope calculation shows that for a Keplerian disk, i.e, k — Q — (GMgtar/?'^)^''^, 
Q ~ 7'-3/2-?/2+p^ where T ~ r"'^ and E ~ represent the power laws for the effec- 
tive temperature and surface density profiles, respectively The typical disk effective 
temperature profile falls off as roughly r~^/^, with a large amount of scatter (Beck- 
with et al. 1990; Kitamura et al. 2002), and disk surface density profiles are beheved 
to lie between about p 1 to 1.5 (Beckwith et al. 1990; DuUemond et al. 2007). 
With q = 1/2 and p = 3/2, one expects Q ~ r^^^^, and at least the outer disk will 
be unstable against GIs. What about GIs in the planet-formation region of the disk, 
presumably r < 40 AU? 

Over the past several decades, the typical picture of a T Tauri disk has been 
the Minimum Mass Solar Nebula model (WeidenschiUing 1977; Hayashi et al. 1985), 
which is based on the known mass distribution of solids in the Solar System and 
which represents the minimum mass, ~ 10~^ Mq, required to form these bodies. 
Beckwith et al. (1990) found that for a sample of 86 T Tauri stars in the Taurus- 
Auriga star forming region, 42% of the systems have detectable protoplanetary disks, 
the disk masses range between 0.001 M© and about 1 M©, and the average disk 
mass is about 0.02 Mq, with massive disks being rare. Andrews & Williams (2005) 
also studied the Taurus-Auriga star formation region and found that 61% of the 153 
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targets have detectable disks, the median disk mass is ~ 0.005 Mq, and less than 
a few percent of disk masses are likely to be gravitationally unstable based on their 
mass estimates and model assumptions. Eisner & Carpenter (2006) found a similar 
median mass for a sample of 336 disks in the Orion Nebular Cluster, and they found 
a fraction of high mass disks consistent with the fraction in Taurus- Auriga. 

With these disk masses, one might expect GIs to play a minor role in T Tauri 
disk evolution, except for possibly r > 100 AU, and expect the MMSN to be a 
fairly accurate description of T Tauri disks. However, there are several potential 
problems with these observations. Grain growth beyond 1 mm will likely result in 
underestimating disk masses because the surface area of the emitting dust grains will 
have changed (Andrews & WiUiams 2005; Hartmann et al. 2006). In addition, as 
argued by Hartmann et al. (2006), inferred mass fluxes from excess UV and blue- 
optical continuum measurements suggest that mass accretion rates are incompatible 
with disk masses and lifetimes, e.g., 0.005 M© with M ~ 10~^ Mq yr"-*^ in a 1-2 Myr 
old disk. With the uncertainty in measured accretion rates, it is unclear how much 
the Md or M measurements are off, but it does indicate that the measurements may 
only be accurate to within a factor of four (Andrews & Williams 2005). Finally, FU 
Orionis events (see below) accrete 0.01 Mq onto the star during their hundred-year- 
long outbursts. Although these disks are somewhere between Class I and II objects, 
one should keep in mind that any one outburst accretes a MMSN, and a typical T 
Tauri disk may go through approximately 10 FU Ori-hke events in less than about 1 
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Myr (Hartmann & Kenyon 1996). 

As suggested by Eisner & Carpenter (2006), the average low disk masses that 
are observed may be more indicative of fast disk evolution than T Tauri disks never 
having a high-mass disk phase. Based on massive disk estimates of several clusters at 
different ages, Eisner & Carpenter argue that disk masses may change by a factor of 
about a few between 0.3 and 2 Myr. In addition, Andrews & WiUiams (2005) argue 
that the ubiquity of planets, where 10% of low-mass stars have detectable planets, 
may indicate disk measurements are systematically underestimated. The majority 
of T Tauri disks with ages of 1 to 2 Myr may very well be too low mass for GIs to 
activate, but GIs should not be discounted during the first Myr of evolution for most 
T Tauri objects and even in the late stages of disk evolution around the high mass T 
Tauri stars (about 1 to 2.5 Mq). 

1.1 The a Disk 

Shakura & Sunyaev (1973) made the ansatz that the turbulent viscosity in an 
accretion disk responsible for carrying angular momentum outward can be approxi- 
mated by 

u — aCgh, (1-2) 

where h is the disk vertical scale height and a < 1 is a parameter that sets the magni- 
tude of the turbulent viscosity. Shakura & Sunyaev did not identify the source of the 
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turbulence, but their heuristic approach allows for analytic solutions to the accretion 
disk problem and relatively easy numerical modeling. In a razor thin Keplerian disk, 
mass transport is described by (e.g., Hartmann 1998; Balbus & Papaloizou 1999) 



In an isothermal Keplerian disk with negligible self-gravity, one can show that p{z) — 
Po exp(— 0^^^/2c^), where q is the isothermal sound speed. By defining the vertical 
scale height h — S/2po, one finds that to within a factor of order unity, h ~ Cj/Q. 
Varying the gas equation of state does not change the factor of order unity signif- 
icantly, and so the appropriate sound speed can be used in place of q. Using this 
relation, equation (1.3) can be rewritten according to a disk theory: 



and so by specifying a, mass transport can be completely described for a given initial 
surface density profile. In a steady state accretion disk, i/E = constant (Hartmann 
1998), and equations (1.3) and (1.4) can be reduced to 




(1.3) 




(1.4) 



M = 37ri/E 37rQ;c^EQ-^ 



(1.5) 
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Note that equation (1.4) indicates that an a disk is an inherently local description 
for mass transport. Because the energy dissipation per unit area in an accretion disk 
can be described by (Pringle 1981) 



GMM , , 

ft « (1.6) 



energy dissipation in an a disk is local as well (see Balbus & Papaloizou 1999). If 
long-range energy transport and/or angular momentum transport are negligible, then 
an OL prescription for disk evolution would be fairly accurate and advantageous, be- 
cause local simulations, e.g., shearing sheets (Balbus & Hawley 1998), would capture 
disk evolution well. Otherwise, the disk behavior will be dissimilar to what a disk 
theory predicts due to long-range angular momentum and energy fluxes (Balbus & 
Papaloizou 1999). Identifying the principal angular momentum transport mechanism 
and how it behaves is critical to understanding whether disks can be described by an 
a disk model. 



1.2 GIs, MRI, and Dead Zones 

There are two mechanisms that have been demonstrated to work efficiently at 
transporting mass inward and angular momentum outward: GIs and the magnetoro- 
tational instability (MRI; see Balbus & Hawley 1991; Desch 2004). As discussed 
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above, GIs require a cold, massive environment to activate. In contrast, the MRI 
in principle only needs a weak magnetic field coupled to ionized species in the gas. 
To see how the MRI produces angular momentum transfer, consider two parcels of 
gas at slightly different radii through which a magnetic field is threaded. As the gas 
orbits, the parcels of gas are separated due to the shear in the disk. Assuming that 
the magnetic fiux cannot diffuse, i.e., d'B/dt = V x v x B (the flux freezing approx- 
imation), the magnetic field remains entrained with the gas. As the gas shears, the 
magnetic fields act as a spring between the two gas parcels. This produces tension 
and an exchange of angular momentum. Because the inner parcel leads in a Keplerian 
disk, the angular momentum transfer is outward, and the mass elements drift further 
apart. This leads to a greater force mediated by the magnetic field lines, and even 
more angular momentum is transfered; an instability ensues (Balbus & Hawley 1998). 

In order for the MRI to activate, ionized species must be present in the gas phase. 
Ionization may be from a thermal source, e.g., coUisional ionization of alkalis, or a 
nonthermal source, e.g., cosmic rays, energetic particles from the star, and X-ray 
irradiation. I refer to these nonthermal sources simply as energetic particles (EPs) 
for this discussion. Thermal ionization of alkalis only occurs for T > 1000 K, and so 
only the innermost portion of a T Tauri disk is expected to be thermally ionized. If 
most of a T Tauri disk is MRI active, it must be due to nonthermal sources. Gammie 
(1996) had the insight that because one expects EPs to be attenuated by the gas, 
there may be regions where the MRI is active and other areas where it is absent. In 
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the inner regions of a disk where the column densities are large, MRI may only be 
active in a thin layer, resulting in layered accretion. As one moves outward in the disk 
and the column density decreases, the entire disk can become MRI active. The region 
where MRI is mostly absent, except for a thin layer at high altitude, is called the 
dead zone. This dead zone is of particular interest to GI studies because mass may 
pile up in dead zones due to the sudden drop in accretion rate as the disk transitions 
from a fully active MRI a disk to a thin, layered accretion flow. 

EPs are attenuated with a scale length of about 100 g cm~^ (Stepinski 1992), and 
so even for a MMSN, the disk will hkely exhibit layered accretion (Desch 2004) and 
a dead zone can form. However, a dead zone may not be tranquil due to Reynolds 
stresses from the thin, MRI-active upper layers (Fleming & Stone 2003). Nonetheless, 
even if mass accretion is only reduced and not altogether halted, mass may still pile 
up in the dead zone (Oishi et al. 2007). If enough mass accumulates, then even for 
an otherwise low-mass disk, GIs can activate. Such mass concentrations may play an 
important role in the FU Ori phenomenon. 

1.3 FU Orionis Phenomenon 

The FU Ori phenomenon is characterized by a rapid (1-lOs yr) increase in optical 
brightness of a young T Tauri object, typically by 5 magnitudes. Emission line spectra 
and strong near and mid infrared excess indicate that the event is driven by sudden 
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mass accretion of the order lO^^M© yr~^ from the inner disk onto the star (Hartmann 
& Kenyon 1996). Because FU Ori objects appear to have decay timescales of about 
100 yr, 0.01 Mq can be accreted onto the star. Note that this is the entire mass of 
the MMSN. 

FU Ori objects are very rare (5-10 known objects, with some objects uncertain; 
Green et al. 2006), but when compared with local star formation rates, it is plausible 
for most T Tauri stars to have approximately 10 FU Ori outbursts with a low state 
lasting 10^ to 10^ yr between outbursts (Hartmann & Kenyon 1996). Some of these 
systems are still surrounded by a substantial remnant envelope (e.g., V1057 Cyg), with 
continued infall onto the disks. Vorobyov & Basu (2005, 2006) find that in their 2D 
magnetohydrodynamics simulations with self-gravity, mass accretion onto a disk from 
an envelope results in episodic bursts of GI activity with about the same frequency 
as expected for FU Ori objects. They suggest that the FU Ori phenomenon may 
be driven by fragmented clumps accreting onto the star. However, Zhu et al. (2007) 
find that this scenario may provide too much disk emission for r > 10 AU to explain 
the FU Ori observations. Moreover, it is unclear whether all FU Ori objects have 
significant envelopes (e.g., FU Ori; Herbig 1977; Kenyon & Hartmann 1991; Green 
et al. 2006; Zhu et al. 2007); the FU Ori phenomenon may not require an infaUing 
envelope to activate. Another external source for driving the FU Ori phenomenon is 
close encounters with a companion. Given the assumption that most T Tauri stars go 
through an FU Ori phase, this appears to be a reasonable idea. Even though several 
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FU Orionis objects are confirmed binaries (e.g., L1551, RNO IB and IC), some show- 
no signs of a nearby companion as determined from tlie absence of spectral line drift 
(Hartmann & Kenyon 1996; but see also arguments by Reipurth 2005 in support of 
a binary-driven mechanism) . 

To date, the best explanation for the optical outburst is a thermal instability. 
Models of this mechanism produce timescales and observational features that are 
similar to the objects that are observed (e.g.. Bell & Lin 1994). For a simple a disk 
model, M Stti/H ~ c^E/Q, and so an increase in the sound speed increases mass 
accretion. If dissipation from mass accretion begins to occur faster than radiative 
cooling, the disk could heat up until hydrogen thermally ionizes. This ionization 
creates an H~ opacity front that strongly reduces the efficiency of radiative cooling. 
The disk continues to heat, and a thermal runaway begins until the unstable region is 
completely ionized. The MRI would likely be active in such a hot disk due to thermal 
ionization, which can lead to a disk-hke accretion (cf. Fromang et al. 2004). Despite 
the success of the thermal instability in explaining the FU Ori phenomenon, there is 
still an unanswered question: What caused the dissipation to become greater than 
radiative cooling? 

Armitage et al. (2001) suggested that GIs in a bursting dead zone might be able 
to trigger an FU Ori outburst by rapidly increasing the accretion into the inner disk 
and initiating a thermal MRI. Likewise, Hartmann (2007, private communication) 
and Zhu et al. (2007) suggest that gravitational torques might drive the temperature 
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near 1 AU above 1000 K and thermally ionize alkalis. The temperatures may need to 
be greater than 1400 K to prevent depletion of ions by dust grains (Sano et al. 2000; 
Desch 2004), but this does not change the general picture. Once the alkalis are ionized, 
a thermal MRI could operate and feed mass inside 0.1 AU until a thermal instability 
activates. The FU Ori phenomenon may be a result of a cascade of instabilities, 
starting with a burst of GI activity in a dead zone, followed by accretion due to a 
thermal MRI, followed finally by a thermal instability. Indeed, recent observations of 
FU Ori indicate that very large mass fluxes are present out to at least r = 0.5 AU 
(Zhu et al. 2007). 

1.4 Chondrules 

Chondrules are small, 0.1-1.0 mm in diameter, igneous globules that were flash 
melted during their formation. They are a fundamental primitive solid inasmuch 
as they can account for over 80% of some chondritic meteorite masses (Hewins et 
al. 2005). Although details are still being debated, most chondrules formed in the 
first 1 to 3 Myr of the Solar Nebula's evolution (Bizzarro et al. 2004; Russell et 
al. 2005). Furthermore, about half of all material accreted onto Earth each year is 
in the form of chondritic meteorites (Hewins et al. 1996), and therefore chondrules 
are arguably the building blocks of the terrestrial worlds. A theoretical description 
of the environments in which chondrules can form is a crucial step in understanding 
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the origin of the Solar System. 

Meteoritics combines astrophysics with geology, and unlike many areas of as- 
trophysics, specimens almost literally fall into our hands. As a result, formation 
constraints for chondrules can be tested through laboratory experiments. Chon- 
drule precursors were flash melted from solidus to liquidus, where high temperatures 
T ~ 1700 K were experienced by the precursors for a few minutes. The melts then 
cooled over hours, with the actual cooling time depending on chondrule type (Scott 
& Krot 2005). Chondrules have diverse petrologies, varying in volatile abundances, 
mineral composition, oxygen isotope ratios, and textures (Jones et al. 2005). Many 
chondrules have flne-grained and igneous coarse-grained rims, and many chondrules 
indicate multiple collisions, as inferred from chondrule fragments and compound chon- 
drules, i.e., chondrules inside chondrules (Scott & Krot 2005). In order to remain 
liquid and preserve volatiles, chondrules are believed to have formed in regions of 
high pressure (10~^ to 10"^ bar) or in a dusty environment, with a dust to gas ratio 
10 to 100 times greater than the typically 1/100 value assumed for the Solar Nebula 
(Wood 1963; Scott & Krot 2005; Hewins et al. 2005). 

Ca-Al-rich Inclusions (CAIs) and Amoeboid Olivine Aggregates (AO As) (e.g., 
MacPherson et al. 2005) also provide constraints on the environment of the early Solar 
Nebula. These refractory inclusions are thought to be among the first solids processed 
in the Solar Nebula, and the formation of CAIs is often used to set the age t = 0. CAI 
and possibly AOA formation occurred in less than 0.3 Myr, and chondrule formation 
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began 0-2 Myr later and lasted for < 3 Myr for normal chondrules (Amelin et al. 2002; 
Itoh et al. 2002; Bizzarro et al. 2004; see Scott & Krott 2005 for a review). CAIs 
formed at pressures near 10~^ atm and at temperatures T ~ 1400 K. Some CAIs were 
completely melted and others were not (MacPherson et al. 2005). Moreover, some 
CAIs show signs of reprocessing, and so they likely experienced chondrule-forming 
events. I do not dismiss the importance of CAI and AOA formation in understanding 
the early conditions of the Solar Nebula. Indeed, they are a critical component, but 
for this dissertation, I focus on chondrule formation. 

Chondritic parent bodies show a range in composition, and this variation has been 
interpreted to be due to temporal and spatial formation differences (Wood 2005). 
There are three basic types of chondrites in which almost all chondrules can be classi- 
fied: carbonaceous, ordinary, and enstatite chondrites. For an in-depth discussion of 
chondrite classification and petrology, I refer the reader to the recent review articles 
by Hewins et al. (2005), Jones et al. (2005), and Scott & Krot (2005). 

Chondrules are separated in their parent bodies by a material called the matrix. 
This material consists of fine-grained minerals and amorphous material, and may be 
related to, if not the same as, the material in fine-grained rims surrounding some chon- 
drules. Generally, chondrules and the matrix show deviations from solar abundance, 
but together are closer to solar values than any component alone, with variation in 
fractionation between chondrule types (Cuzzi et al. 2005; Huss et al. 2005). This 
abundance correspondence is referred to as chondrule-matrix complementarity. 



1.4. CHONDRULES 



17 



Chondrule coUisional histories, isotopic fractionation data (see below), chondrule- 
matrix complementarity, fine-grained rim accumulation (Morfill et al. 1998), and 
petrological and parent body location arguments (e.g.. Wood 2005) suggest that chon- 
drules formed in the Solar Nebula, as suggested by Wood (1963), in strong, localized, 
repeatable heating events. For the discussions in this dissertation, I will focus on 
the shock wave model for chondrule formation, because it is the most well-developed 
and viable formation hypothesis (lida et al. 2001; Desch & Connolly 2002; Cielsa 
& Hood 2002; Boss & Durisen 2005a,b; Hartmann 2005; Miura & Nakamoto 2006). 
However, other formation mechanisms may play a role in processing some solids, e.g., 
the X-Wind model (Shu et al. 2001) and hghtning (Whipple 1966; Pihpp et al. 1998; 
Desch & Cuzzi 2000). It should also be noted that the shock wave model does not 
necessarily identify the shock-driving mechanism, and identifying this mechanism is 
an ongoing topic for investigation. 

Shock waves can form chondrules through gas-drag heating. The stopping time 
for a particle with radius a < 1 cm 



Pad 



(1.7) 




where pa is the particle density and p the gas density (Cuzzi et al. 2001). As the 
stopping time becomes very short, it is possible to melt chondrule precursors for a 
range of pre-shock conditions through a combination of friction and interaction with 
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the hot, post-shock gas (Desch & Connolly 2002; Ciesla & Hood 2002). The melts' 
interaction with each other and the gas prevents the newly-forming chondrules from 
cooling too quickly 

A critical chondrule-formation constraint is that chondrules are depleted in some 
elements, e.g., K, Fe, Si, Mg, presumably from devolatization during the precusor 
melting, but the corresponding isotopic fractionation of species in chondrules, e.g, 
sulfur in troilite, is mostly absent (Tachibana & Huss 2005). There are two favored 
explanations for this problem: (1) Chondrule heating is so rapid that isotopic frac- 
tionation is suppressed. (2) Chondrules were embedded within the evaporated gas 
of other chondrules, and came to equilibrium with that gas. Miura & Nakamoto 
(2006) demonstrate through ID radiation-shock calculations that the optical depth 
as measured between the shock front to about 10^° cm upstream cannot exceed r ~ 
1-10, depending on the density and the pre-shock velocity. They found that unless 
this criterion is met, chondrule precursors spend too much time above T ~ 1400 K to 
avoid isotopic fractionation of sulfur. However, laboratory experiments indicate that 
the rapid heating and cooling needed to avoid isotopic fractionation is incompatible 
with observed chondrule textures (Hewins et al. 2005), and so the criterion may not 
be valid. In contrast, Cuzzi & Alexander (2006) argue that chondrules can equili- 
brate with their surroundings if the chondrule-forming region is larger than a few 
thousand kilometers because the evaporated gas cannot diffuse in the time it takes 
for the chondrule to be processed. 
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Another point to consider is that chondrule formation may not be hmited to 
asteroid belt distances. Stardust results (e.g, McKeegan et al. 2006) as well as comet 
observations (Wooden et al. 2005) indicate that high-temperature thermally processed 
solids are present in planetesimals from the outer Solar System (comets) , and so either 
large-scale radial transport occurred in the Solar Nebula or some high-temperature 
processing occurred in the outer disk. Understanding how chondrules form and how 
processed solids are mixed in the Solar Nebula is critical to describing its evolution, 
and ultimately, planet formation. 

One plausible shock-driving mechanism is a global spiral wave (Wood 1996). 
Harker & Desch (2002) suggest that spiral waves could also explain thermal pro- 
cessing at distances as large as 10 AU, and Boss & Durisen (2005a,b) suggest that 
GIs may be able to produce the required shock strengths to form chondrules. In 
addition, Boley et al. (2005) suggest that spiral waves could be a source of turbulence 
as well as a shock mechanism. Global spiral shocks are appealing because they fit 
many of the constraints above. They may be repeatable, depending on the formation 
mechanism for the spiral waves; they are global but produce fairly local heating; they 
can form chondrules in the disk; and they can work in the inner disk as well as the 
outer disk. 
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1.5 Disk Fragmentation and the Planet Formation 
Debate 

Knowing under what conditions protoplanetary disks can fragment is crucial to 
understanding disk evolution inasmuch as a fragmented disk may produce gravita- 
tionally bound clumps. This is the disk instability hypothesis for the formation of 
gas giant planets (Kuiper 1951; Cameron 1978; Boss 1997, 1998). The strength of 
GIs is regulated by the cooling rate in disks (Tomley et al. 1991, 1994; Pickett et 
al. 1998, 2000a, 2003), and if the cooling rate is high enough in a low-Q disk, a disk 
can fragment (Gammie 2001). Consider the two-dimensional adiabatic index F, such 
that J pdz — P — XE'", where K is the poly tropic constant. Gammie quantified 
that for r = 2, a disk will fragment when tcooi^ ^ 3. Here, icooi is the local cooling 
time and Q is the angular speed of the gas. This criterion was approximately con- 
firmed in 3D disk simulations by Rice et al. (2003) and Mejia et al. (2005). Rice et 
al. (2005) showed through 3D disk simulations that this fragmentation criterion de- 
pends on the adiabatic index and, for 7 = 5/3 or 7/5, the fragmentation limit occurs 
when tcooi^ ^ 6 or 12, respectively. These results show that a change by a factor 
of about 1.2 in 7 has a factor of two effect on the critical coohng time. In addition, 
these results indicate that the cooling time must be roughly equal to the dynamical 
time of the gas for the disk to be unstable against fragmentation when 7 = 5/3. Do 
these prodigious cooling rates occur in disks when realistic opacities are used with 
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self-consistent radiation physics? 

Nelson et al. (2000) used 2D SPH simulations with radiation physics to study 
protoplanetary disk evolution. Because their simulations were evolved in 2D, they 
assumed that the disk at any given moment was in vertical hydrostatic equilibrium. 
Using a polytropic vertical density structure and Pollack et al. (1994) opacities, they 
cooled each particle according to an appropriate effective temperature. In their sim- 
ulations, the cooling rates are too low for fragmentation. In contrast. Boss (2001, 
2005) employed radiative diffusion in his 3D grid-based code, and fragmentation oc- 
curs in his simulated disks. Besides the difference in dimensionality of the simulations. 
Boss assumed a fixed temperature structure for Rosseland mean optical depths less 
than 10, as measured along the radial coordinate. Boss (2002) found that the frag- 
mentation in his disks is insensitive to the metallicity of the gas and attributed this 
independence to fast coohng by convection (Boss 2004a). However, it must be noted 
that Nelson et al. (2000) assumed a vertically polytropic density structure. Because 
the entropy S ^ InK, where p = Kp^ is the polytropic equation of state, the Nelson 
et al. approximation assumes efficient convection. One of the principal objectives in 
this dissertation is to understand and resolve these discrepant results. 
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1.6 About This Dissertation 

In this dissertation, I explore the three-dimensional behavior of spiral waves in 
gravitationally unstable disks, and discuss their effects on vertical disk structure, 
mass transport, the FU Orionis phenomenon, and chondrule formation. In addition, 
I investigate under what conditions, if any, unstable disks fragment, whether sudden 
vertical motions in these unstable disks can be better explained by shock bores than by 
convection, whether the details of the treatment for radiation transport are important 
to disk evolution, and whether dust settling affects the strength of spiral shocks. 
The culmination of this work is an effort to understand to what degree chondrule 
formation, planet formation, the FU Ori phenomenon, dead zones, and GI activity 
are interdependent. 

In Chapter 2, I discuss numerical methods that are employed to assess shock 
strengths in these disks, to approximate radiative cooling, and to account for the 
rotational and vibrational states of molecular hydrogen. I also outline the hydrody- 
namics codes used for these studies. In order to test the radiation physics described in 
Chapter 2, several radiation hydrodynamics tests along with their results for the radi- 
ation algorithms implemented here are discussed in Chapter 3. The initial models for 
these studies are detailed in Chapter 4, and in Chapter 5, I outline shock bore theory. 
The role of convection, the locality of mass transport, the importance of radiation 
boundary conditions, and radiative cooling rates in Gl-active disks are examined in 
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Chapter 6. In Chapter 7, I explore the connection between massive, bursting dead 
zones, the FU Ori phenomenon, and chondrule formation. Finally in Chapter 8, I 
summarize the principal conclusions for the work presented in this dissertation. 



Chapter 2 



NUMERICAL METHODS 

There are two versions of the lU hydrodynamics code that have been used for the 
studies described in this dissertation. I outhne both versions in this Chapter, and 
layout algorithms and analysis routines that I have developed for these studies. 

2.1 The Standard lU Code 

The standard version of the lU code (SV) is an explicit, Eulerian code that solves 
the equations of hydrodynamics in conservative form on an evenly spaced, cylindrical 
grid. The code is second-order accurate in space and in time, and it solves for self- 
gravity after the first source half-step (see below). The SV is based on the Tohhne 
(1980) code. Yang (1992) made the code second order in space by including a van Leer 
advection scheme (van Albada et al. 1982), made the code second order in time, and 
modified the mesh staggering (see below) to its current form. Pickett (1995) included 
an energy equation, artificial viscosity (AV), and ad hoc cooling prescriptions. Mejfa 

(2004) added radiation physics, and Cai (2006) made several improvements to the 
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Mejia radiation algorithm. 

The equations of hydrodynamics with self-gravity are 
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dp 
di 



+ V-pv = 0, 



(2.1) 



+ V • pvv = -Vp - pV$ - V • (pQ) , 



(2.2) 



dt 



de 
di 



+ V • ve = -pV • V + r + A, and 



(2.3) 



(2.4) 



i.e., the mass continuity equation, the equation of motion, the energy equation, and 
Poisson's equation, respectively. Here, p is the mass density, v is the velocity of the 
fluid, p is the gas pressure, $ is the total gravitational potential, and e = pe, where 
e is the specific internal energy of the gas. The last term in the equation of motion 
is the artificial viscosity term. The tensor Qij = CQ{Avi)'^ for j — i and Avi < 0, 
but otherwise, as in the von Neumann & Richtmeyer scheme (Norman & Winkler 
1986). The constant Cq is an experimentally determined parameter that spreads 
shocks over several cells to include the appropriate amount of entropy generation and 
to stabilize shocks against ringing (see Pickett 1995). The AV contribution to the 
energy equation is given by P = p {QrrdrVr + Q(j)(j)r~^dri,Vcj) + QzzdzVz), and radiation 
cooling and heating are described by the A term (see §2.4). If an equation of state is 
chosen such that p = f (p), e.g., barotropic gases, only the mass continuity equation 
and the equation of motion are required to describe the flow. However, \ip — f (p, e). 
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the energy equation is also required. The SV version of the code employs an ideal 
gas equation of state with a constant ratio of specific heats 7 such that 



p = (7 - 1) e. (2.5) 

The simple form for p in equation (2.5) allows the energy equation to be recast in 
conservative form (Williams 1988), namely 

^ + V • ve^T = (r + A) . (2.6) 

at 7 

Pickett (1995) implemented equation (2.6) in the SV and found that it produced 
slightly better results than using equation (2.3). 

In addition to the mass density and the internal energy density of the fluid, the 
SV uses the radial momentum density S, the vertical momentum density T, and the 
angular momentum density A to evolve the flow. The array values for p, e, and A are 
set at cell centers, while the S and T array values are set on the center of cell faces 
(see Fig. 2.1). Yang (1992) found that this staggered grid formalism is more stable 
than a complete cell-centered formalism. 

Because a cylindrical grid is used, equations (2.1-2.4) are solved in cylindrical 
coordinates, and so equation (2.2), in particular, needs to be rewritten in a more 
useful form. Let pv = Sir + A/re^ + Tiz. The divergence term in equation (2.2) can 



2.1. THE STANDARD lU CODE 



27 



be written as 

V • pvv = pvV • V + V • Vpv, (2.7) 
for position vector r = re^ + ze^- The radial component of equation (2.7) is 

(V • pvv) • er = r'^dr (rSvr) + r'^d^ (Sv^) + (Sv^) - r'^v^A 

= V-Sv-^. (2.8) 

Likewise, the azimuthal component can be written 

r (V • pvv) • = r~^dr (rAvr) + r~^d^ i^v^) + (Avz) 

= V-Av. (2.9) 

The vertical component is straightforward, and the equation of motion can now be 
written as 



^ + V5v = - ^- drpQrr ^ A' 

dt dr dr rdr pr^^ 

dT ^ ^ dp a$ dpQ,, 

_ + v.rv = -,£-,--J^. (2.12) 
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Figure 2.1 Location of values for the five principal hydrodynamics arrays. The index- 
ing for the cyhndrical grid (r, z, 0) = (J,K,L). 
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The SV uses an operator splitting method to evolve the hydrodynamics. The two 
operations are called sourcing, S, which advances the right hand side in equations 
(2.10-2.12), and fluxing, JF, which advances the advection terms in equations (2.1), 
(2.3), and (2.10-2.12). The solution for a hydrodynamic variable X for any one time 
step is dX/ dt = S+T. As described in Yang (1992), the fluxing is determined through 
a van Leer advection scheme (van Albada et al. 1982), with the fluxed quantities being 
S, T, A, p, and e^T. For sourcing the arrays, three separate calculations are made. 
First, the forces due to pressure gradients and potential gradients are determined by 
center-differencing schemes and used to update the momentum densities. Second, 
the effects of AV on the momentum densities are included. Finally, heating due to 
AV and cooling due to either ad hoc prescriptions or radiation transport are used to 
update the internal energy densities. It is possible to run the simulation with cold AV 
(Pickett & Durisen 2007), i.e., AV only affects the equation of motion. However, this 
is only done under special circumstances. To illustrate how the fluxing and sourcing is 
used to achieve second-order time integration, I reproduce the flow chart from Mejia 
(2004, Fig. 2.4) in Figure 2.2. Note that the pressure is updated before each sourcing 
and the potential is calculated only once before the second sourcing. 



30 



CHAPTER 2. NUMERICAL METHODS 




Figure 2.2 A flow chart illustrating the second-order time integration scheme. Repro- 
duced from Mejia (2004), Figure 2.4. 
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The total potential $ = + where $6 is some background potential and $p 
is the potential due to the mass on the grid. The background potential can be set 
to zero, as is done for protostar+disk models (e.g., Yang 1992; Pickett 1995; Pickett 
et al. 1996, 1998). For the models presented here and for a multitude of previous 
studies (Pickett et al. 2001, 2003; Mejia 2004; Mejia et al. 2005; Cai et al. 2006, 2007; 
Pickett & Durisen 2007), only the disk is evolved, and so the background potential 
is a point mass or fixed potential due to the star that is otherwise excluded from 
the simulations. This background potential is held fixed at the center of the grid. 
Although keeping the star rigid may create some spurious one-arm behavior in the 
simulations, simply moving the star to locations that stringently keep the center of 
mass at the center of the grid (see Boss 1998) may be problematic as well. Routines 
are being developed to explicitly model the motion of the star, but have not been 
employed for the studies discussed here. 

The potential due to the mass is found in two parts. First, the potential on the grid 
boundary is calculated through a spherical harmonics expansion with / = \m\ = 10. 
The boundary solution then provides the Dirichlet boundary conditions for a direct 
Poisson solver (Tohline 1980). The Poisson solver uses a discrete Fourier transform 
to decompose the 3D potential problem into LMAX 2D problems, for grid dimensions 
(r, z, (p) = (JMAX, KMAX, LMAX). The series of 2D problems can then be solved in 
parallel with a cyclic reduction method (Tohline 1980; Pickett et al. 1996). Pickett et 
al. (2003) demonstrated that when a constant density blob was loaded onto the grid. 
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the Z = |m| = 10 expansion is sufficient to describe tfie boundary potential and tfiat 
tlie solution's accuracy is determined by grid resolution, not the boundary potential 
expansion. Moreover, I have used the Cohl & Tohline (1999) Bessel expansion to 
test the accuracy of the SV boundary potential solver for the Wengen test 4, which 
is a set of disk simulations by a wide variety of codes for the same initial conditions 
(see Mayer et al. 2007, in preparation). I find that the outcome is unchanged when 
the Legendre half-integer polynomials are carried out to m = 30, where the Bessel 
expansion showed convergence to 1 part in 10 million when m was increased to 40. 
The boundary potential solution reaches the required accuracy with I — |m| = 10. I 
note that the Wengen test 4 is challenging for a potential solver because the disk is 
very flat, asymmetric, and develops dense clumps. The potential solver is sufficient 
for these tests. 

2.2 CHYMERA 

The hydrodynamics routines in CHYMERA, Computational HYdrodynamics with 
MultiplE Radiation Algorithms, are different from those in the SV in one significant 
way: the rotational and vibrational states of molecular hydrogen are explicitly cal- 
culated (see §2.3), and so a 7 = constant approximation is no longer valid. This, 
in turn, means that the energy equation cannot be recast into its conservative form 
(equation [2.3]), and so e must be fluxed and — pV ■ v must be explicitly calculated 
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during sourcing. For the sourcing of the energy equation, the work term is calculated 



where only the relevant indices are shown, all velocities are face-centered, and p is 
the average of the current half step's pressure and a provisional pressure (Black & 
Bodenheimer 1975), i.e., the pressure if the internal energy were updated according to 
p — p. Without the inclusion of the provisional pressure, CHYMERA cannot match 
the Sod shock tube (Sod 1978; Hawley et al. 1984) results of the SV when the shock is 
forced to propagate in the vertical direction. However, when the provisional pressure 
is included, the solutions agree to within a few hundredths of a percent, and energy 
loss is kept to about a tenth of a percent (Figure 2.3; see Pickett 1995, Appendix B 
for test details) . 

In addition to changing the energy equation, the equation of state must also be 
altered. Instead oi p = (7 — l)e, p = kpT/fj,mp, where T is interpolated from a 
table based on a given cell's e and where /i is the mean molecular weight specified for 
calculating the e — T table. 



by 




+ 



t'^(L + 1) - Vct>{L) 



(2.13) 



+ 
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Figure 2.3 Results for the Sod shock tube test for the SV (squares) and CHYMERA 
(Xs). The setup is two uniform regions of 7 = 7/5 gas in contact with each other and 
initially at rest. One region is a p = 1 g cm^'^, p = 1 dyne cm^^ gas, and the other 
a p — 0.125 g cm~^, p = 0.1 dyne cm"^ gas. The plots for density, velocity, pressure, 
and specific energy, in cgs units, are compared with analytic values (solid curves). 
The results are, for the most part, indistinguishable when the shock propagates along 
the vertical direciton. 
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2.3 H2 Thermodynamics 

As discussed in the Introduction, the disparities between disk evolutions reported 
by various hydrodynamics groups are hkely due to a combination between different 
treatments of radiation transfer and the rotational states of molecular hydrogen. Be- 
cause the thermodynamics can strongly change the outcome of a simulation (Pickett 
et al. 1998, 2000), I have implemented in CHYMERA an internal energy that takes 
into account the translational, rotational, and vibrational states of H2. During its 
development, Boley et al. (2007a) noticed that in all treatments to date of planet 
formation by disk instability, the effects of the rotational states of H2 have been, at 
best, only poorly approximated, and they drew attention to possible consequences of 
various approximations for the internal energy of H2 that are in the literature. I out- 
hne the algorithm employed in CHYMERA, and recapitulate several of the cautions 
presented by Boley et al. (2007a). 

For this discussion, I refer the reader to Pathria (1996). Consider the following 
thermodynamic properties of an ideal gas: Let E be the internal energy for N par- 
ticles, e the specific internal energy, e the internal energy density, p the pressure, 
T the gas temperature, p the gas density, /x the mean molecular weight in proton 
masses, Cy the specific heat capacity at constant volume, Z the partition function for 
the ensemble, z the partition function for a single particle, and R = k/nip, where 
k is Boltzmann's constant and nip is the proton mass. I only consider independent 
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contributions to the partition function from translation, rotation, and vibration rep- 
resented by Z = = = (^tran2;rot^vib)^- The internal energy E and the 
specific internal energy e can be calculated by 



E ^ NkT^——- e= -T^—— (2.14) 



for constant p. Because the gas is ideal, = de/dT. 

2.3.1 Molecular Hydrogen 

Molecular hydrogen exists as parahydrogen and as orthohydrogen where the pro- 
ton spins are antiparallel and parallel, respectively. The partition function for parahy- 
drogen is 

Zp^Yl + 1) exp (-J (j + 1) e^jT) , (2.15) 

Jeven 

and the partition function for orthohydrogen is 

Zo = J2^ (2j + 1) exp i-j (j + 1) e^jT) , (2.16) 

where ^rot — 85.4 K (Black & Bodenheimer 1975). When the two species are in 
equilibrium, ^rot — Zp + Zo- However, the ortho/para ratio (b:a) could also be frozen 
if no efficient mechanism for converting between the species is available. This leads 
to Zrot — Zp"'^^"''^''^^ z'o''^^'^~^''^\ where z'^ — Zgexp (26,:ot/T). The additional exponential 
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is required in the orthohydrogen partition function when the ortho and para species 
are at some fixed ratio to ensure that rotation only contributes to the internal energy 
once the rotational states are excited, i.e., z'^ — > constant as T — > 0. 

To consider the vibrational states, I approximate the molecule as an infinitely 
deep harmonic oscillator, where 



^vib = 5987 K (Draine et al. 1983). At this time, CHYMERA is only designed 
to investigate temperatures T < 1500 K where dissociation of H2 is insignificant. 
For this reason, 1 choose to ignore differences between equation (2.17) and a proper 
Zyih, which would take into account the anharmonicity of the molecule and that the 
molecule has a finite number of vibrationally excited states. 

I can use equation (2.14) to write the specific internal energy for H2 



The specific internal energy profiles as calculated by equation (2.18) are shown in 
Figure 2.4 (left panel) for an equilibrium mixture (solid black), pure parahydrogen 
(solid red), and a 3:1 ortho/para mixture (solid blue). The offset of the correct 3:1 
mix profile is due to the energy stored in the parallel spins of the protons. When 



1 



(2.17) 



l-exp(-^vib/T) 




(2.18) 
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the gas is ideal and dissociation and ionization can be ignored, the first adiabatic 
exponent Fi = 1 + R/nCy = Cp/cy — 7 (Cox & Giuh 1968). Sohd curves in Figure 
2.4 (right panel) indicate the Fi profiles for the corresponding solid curves in the left 
panel, which is consistent with Figure 2 of Decampli et al. (1978), as it should be 
because my derivation of is equivalent to theirs. All dashed curves and the curves 
in the center panel result from various approximations for e, as discussed below. 

2.3.2 Approximations for e 

There are several approximations for e that are employed in the literature, and 
several of these approximations can have strong consequences for disk dynamics, as 
argued by Boley et al. (2007a). If is constant, then e = CyT. Boley et al. (2007a) 
pointed out that Black & Bodenheimer (1975) calculate c„ from the Helmholtz free 
energy, which is valid, but then assume e = c^T, which is invalid, because is 
dependent on T. Other authors have followed suit (e.g., Whitehouse & Bate 2006; 
Stamatellos et al. 2007), and this assumption could be troublesome for gas dynamics 
in a hydrodynamics simulation. Curves for c„T are also shown in Figure 2.4 (left 
panel) by dashed lines. 
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1 1.5 2 2.5 

log T (K) 

Figure 2.4 Left panel: The solid curves show the proper specific internal energy for 
an equilibrium mixture (black), a pure parahydrogen gas (red), and a constant 3:1 
ortho/parahydrogen mixture (blue). The dashed curves indicate e = c^T, which can 
deviate strongly from the correct curves. Right panel: Fi curves for the corresponding 
solid curves in the left panel. The dashed curve demonstrates the consequences for 
Fi when a seemingly innocent interpolation is used for e (see Boley et al. 2007b). 
Center panel: Corresponding Fi curves for the dashed curves in the left panel. The 
deviations from the correct e have strong consequences for the Fi profiles. 
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The approximation e = CyT gives quite different behavior from the correct e, 
e.g., the incorrect 3:1 curve most closely follows the correct pure parahydrogen curve. 
The dynamical effects that could result from assuming e = CyT are evaluated by 
taking the temperature derivative of the dashed curves in Figure 2.4 (left panel). The 
resulting Fi profiles are shown in the center panel of Figure 2.4, and these curves 
are very different from the Fi profiles that they are meant to approximate (right 
panel). As mentioned above, Cy can be calculated from the Helmholtz free energy 
as is done by Black & Bodenheimer (1975). This makes it possible to compute Fi 
from Cy correctly but then evolve the gas with an erroneous effective specific heat 
because many hydrodynamics codes evolve e = pe (e.g.. Black & Bodenheimer 1975; 
Boss 1984a, 2001; Monaghan 1992; Stone & Norman 1992; Pickett 1995; Wadsley et 
al. 2004). If the ideal gas law p — pRTf/i is assumed as well, then effective Fi profiles 
hke those shown in Figure 2.4 (center panel) seem to be unavoidable when e = CyT is 
assumed for a temperature dependent Cy. Because fragmentation becomes more likely 
as Fi becomes smaller (Rice et al. 2005), the e = CyT assumption should artificially 
make fragmentation more likely in some temperature regimes and less likely in others. 
Boley et al. do note, however, that the severity of this error may depend on the state 
variables evolved in a given code, and only the authors who employ e = CyT will be 
able to say in detail how it affects their simulations. 

Boley et al. (2007a,b) also caution against simple approximations to e because 
they too may have strong consequences for Fi. For example. Boss (2007) uses a 
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quadratic interpolation for e between 100 K and 200 K.^ In Figure 2.4 (right panel, 
dashed curve), I illustrate the consequences for Fi when using this e approximation 
by calculating the specific heat at constant volume c„ = de/dT and by relating that 
to Fi (see Cox & Giuli 1968). Because fragmentation becomes more likely for lower 
Fis (Tomley et al. 1991; Boss 1997, 2000; Pickett 1998; Rice et al. 2005), an approx- 
imation for e like that used by Boss between 100 and 200 K is likely to make the 
disk susceptible to fragmentation for that temperature range. Finally, a constant Fi 
approximation (Pickett et al. 2003; Rice et al. 2003; Lodato & Rice 2004; Mayer et 
al. 2004, 2007; Rice et al. 2005; Mejia et al. 2005; Cai et al. 2006; Boley et al. 2006) 
poorly represents e between about 60 and 300 K, a plausible temperature regime for 
the formation of Jupiter; neither Fi = 5/3 nor 7/5 can be assumed confidently. 

Preliminary simulations of a disk with solar composition indicate that when GIs 
activate between 30 and 50 K for an equilibrium ortho/para mixture, the e = c^T 
simulation evolves more rapidly and has a more flocculent spiral structure than the 
correct e simulation for the same cooling rates. In addition, denser substructures 
form in some spiral arms of the e = CyT simulation throughout the simulation, while 

dense substructures only form during the burst of the correct e simulation^. When 

-'^Bolcy et al. (2007a) stated that Boss uses a discontinuous internal energy for H2, where H2 
contributes SkT/^nip to the specific internal energy e for T < 100 K and ^kT/ArUp for T > 100 K, 
based on his own citations to Boss (1984a). Through private communication, Boss indicated that 
he now uses a quadratic interpolation for e. This change is not mentioned in the literature, but 
according to Boss, the interpolation has been used in all his simulations since Boss (1989). Boley et 
al. (2007b) clarified this discrepancy. 

^The evolution of these simulations can be viewed at http:/hydro. astro. indiana.edu/westworld. 
Click on the link titled "H2 ortho-para equilibrium tests" under the "Movies" tab. 
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the instabilities occur outside this temperature regime, the differences are diminished. 



2.3.3 The Ortho/Para Ratio 

As indicated in the previous section, the dynamical behavior of the gas is depen- 
dent on the ortho/para ratio and whether the species are in equilibrium for all T. 
This ratio for various astrophysical conditions has been addressed by several authors 
(e.g., Osterbrock 1962; Dalgarno et al. 1973; Decamph et al. 1978; Flower & Watt 
1984; Sternberg & Neufeld 1999; Fuente et al. 1999; Rodriguez-Fernandez et al. 2000; 
Flower et al. 2006) , typically in the context of interstellar clouds or photodissociation 
regions. However, for plausible Solar Nebula conditions, the ortho/para ratio has 
been inadequately addressed. For example, Decamph et al. (1978) used an estimate 
for the H+ number density that was derived originally to give the total gas phase ion 
number density in gas for which dissociative recombination dominates the removal of 
ions. At protoplanetary disk number densities, however, ion removal should be pri- 
marily on grain surfaces if the ratio of grain surface area to hydrogen nucleon number 
density is the same as it is in diffuse interstellar clouds. 

For protoplanetary disk conditions, the conversion between ortho and parahydro- 
gen is principally due to protonated ions such as Hg . Consequently, I will assume 
that all ionizations lead to Hg formation. Another possible conversion mechanism is 
through interactions between H2 and grains. However, this conversion might only be 
significant when the temperature drops below about 30 K (Le Bourlot 2000). 
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Consider the balance between production by cosmic rays (CR) and deple- 
tion by dust grains: 

(n {B.2) = 'ngTra^vrii, (2-19) 

where n (H2), n^, and Ug are the H2, H3 , and grain number densities, respectively, ( is 
the ionization rate by CRs and other energetic particles (EP), a is the average radius 
of the grains, and v is the thermal velocity of H3 . If standard interstellar extinction is 
assumed, then a = n (H2) /ngiia^ ~ 10^^ cm~^, but as discussed below, this number 
is ambiguous. The C appropriate for a protoplanetary disk is also ambiguous. Cosmic 
rays and stellar EPs are important in ionizing the disk surface (Desch 2004; DuUemond 
et al. 2007), but because these particles are attenuated exponentially with a scale 
length of about 100 g cm~^ (Umebayashi & Nakano 1981), stellar EPs probably do not 
contribute to nj. Moreover, protostellar winds could lead to a significant reduction of ( 
in analogy to CR modulation by the solar wind (Webber 1998). However, at a surface 
density of roughly 380 g cm~^, EP production by ^^Al decay is as important as CRs, 
with C ~ 10-^*^ s-i (Stepinski 1992). It is likely that 10"^^ s"^ < C < 10"^^ s-^ For 
the following estimate, I adopt the interstellar rate ( — 10~^^ s~^ (Spitzer & Tomasko 
1968). Using these numbers in equation (2.19) and adopting a thermal velocity of 1 
km s~^, rii fa 0.1 cm~^. By adopting a coUisional rate coefficient a = 1 x 10~^ cm^ s"~^ 
for the H3" interaction with H2 (Walmsley et al. 2004), the lower limit timescale for 
ortho and parahydrogen to reach equilibrium is te = {ani)~^ — 300 yr. 
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The equilibrium timescale is short enough that the ortho/para ratio can thermal- 
ize in the hfetime of a disk, but the equihbrium timescale is longer than the dynamical 
timescale inside about 40 AU: ortho and parahydrogen should be treated as indepen- 
dent species for hydrodynamical simulations of young protoplanetary disks. 

What ortho/para ratio should a dynamicist assume for gravitationally unstable 
protoplanetary disk simulations? The answer is uncertain. Vertical and radial stirring 
induced by shock bores (Boley & Durisen 2006, see Chapter 5), which could possibly 
lead to mixing of the low altitude disk interior with the high-altitude photodissocia- 
tion region in the disk atmosphere (DuUemond et a. 2007), will transport gas through 
different temperature regimes on dynamic timescales. This could lead to nonthermal- 
ized ortho/para ratios like those that are measured from H2 rotational transition 
lines in some photodissociation regions (Fuente et al. 1999; Rodriguez-Fernandez et 
al. 2000) and in Neptune's stratosphere (Fouchet et al. 2003). Moreover, accretion of 
the outer disk will bring material with a cold history into warmer regions of the disk. 
It is unclear whether the ortho/para ratio of, say, 15 K gas will be thermalized with 
Zo/ Zp ^ or whether the ortho/para ratio will be 3:1, which is the expected ratio for 
H2 formation on cold grains (Flower et al. 2006). Unfortunately, the ortho/para ratio 
may be critical to the evolution of a protoplanetary disk. As can be seen in Figure 2.4, 
the pure parahydrogen mix has a Fi that approaches 4/3 for T f=:i 160 K. This could 
make the 160 K regime the most likely region of the disk to fragment because, as Fi 
decreases, it becomes harder for the gas to support itself against local gravitational 
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and hydrodynamic stresses (Rice et al. 2005). Hydrodynamicists need to consider 
ortho/para ratios between pure parahydrogen and 3:1 because of our ignorance of 
this ratio in protoplanetary disks. 

The above discussion is based on the assumption that a 10^^ cm~^, which is 
probably reasonable for very young protoplanetary disks but may not be reasonable 
for disks with ages of about 1 Myr. Grain growth and dust settling may significantly 
lower the value of a by depleting the total grain area (e.g., Sano et al. 2000). Because 
models of T Tauri disks must take into account the effects of grain growth in order to 
match observed spectral energy distributions (D'Alessio et al. 2001, 2006; Furlan et 
al. 2006), there may be a period in a disk's evolution when the ortho and parahydrogen 
change from dynamically independent species to species in statistical equilibrium. 
Such a transition may also take place at certain radii in a disk, e.g., near edges of 
a dead zone (Gammie 1996). As indicated by Figure 2.4, a transition to statistical 
equilibrium could have significant dynamical consequences for disk evolution and may 
induce clump formation by GIs. 

2.4 Radiation Algorithms 

2.4.1 The M2004 and C2006 Schemes 

In this section I briefly describe the radiation algorithms developed by Mejia (2004; 
hereafter M2004) and Cai (2006; hereafter C2006), which are also described briefly 
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in Boley et al. (2006, 2007c). In the M2004 scheme, flux-hmited diffusion is used in 
the r, 4>, and z directions on the cyhndrical grid everywhere that the vertically inte- 
grated Rosseland optical depth r > 2/3, which defines the disk's interior. For mass 
at lower optical depths, which defines the disk's atmosphere, the gas is allowed to 
radiate as much as its emissivity allows, with the Planck mean opacity used instead 
of the Rosseland mean opacity. The disk interior and atmosphere are coupled with 
an Eddington-like boundary condition over one cell. This boundary condition defines 
the fiux leaving the interior, which can be partly absorbed by the overlaying atmo- 
sphere. Likewise, feedback from the atmosphere is explicitly used when solving for 
the boundary flux. However, cell-to-cell radiative coupling is not explicitly modeled 
in the disk's atmosphere. This method allows for a self-consistent boundary condition 
that can evolve with the rest of the disk. The C2006 routine (see also Cai et al. 2006) 
improves the stability of the routine, as described in Chapter 3, by extending the 
interior/ atmosphere fit over two cells. 

A problem with the M2004 and C2006 routines (see Chapter 3) is a sudden drop in 
the temperature profile where r = 2/3. The drop is due to the omission of complete 
cell-to-cell coupling in the optically thin regime (r < 2/3). However, as shown in Bo- 
ley et al. (2006; see Appendix B), the boundary does permit the correct fiux through 
the disk's interior. Because the fiux through the disk is correct, the temperature drop 
is mainly a dynamic concern inasmuch as it might seed convection (Boley et al. 2006). 
In order to obtain the correct fiux and temperature profiles, a method for calculating 
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fluxes that takes into account the long-range effects of radiative transfer is required. 



2.4.2 The BDNL Scheme 

In order to account for the long-range effects of radiative transfer, at least in part, 
I have developed a radiation algorithm that couples the flux-limited diffusion used in 
the M2004 and C2006 routines with vertical rays. The method was flrst described 
by Boley et al. (2007c), and so I refer to it as the BDNL scheme. Consider some 
column in a disk with fixed r and 4>. Take that column out of context, and imagine 
that it is part of a plane-parallel atmosphere. In this case, heating and cooling by 
radiation can easily be described with the method of discrete ordinates (see, e.g., 
Chandrasekhar 1960; Mihalas & Weibel-Mihalas 1986). This method uses discrete 
angles that best sample the solid angle, as determined by Gaussian quadrature. In 
a plane-parallel atmosphere, a single ray can provide decent accuracy if the cosine 
of the angle measured downward from the vertical to the ray is fj, — I use 

this approach to approximate radiative transfer in the vertical direction, and include 
flux-hmited diffusion (Bodenheimer et al. 1990) in the r and directions everywhere 
that T > l/y/S. Naturally, this is only a crude approximation when one places the 
column back into context. However, Boley et al. (2007c) argue that this method 
represents the best implementation of radiative physics for simulating protoplanetary 
disks with three-dimensional hydrodynamics thus far, because it captures the long- 
range effects of radiative transfer that are excluded in pure flux-limited diffusion 
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routines and handles optically thick and thin regions. As demonstrated by Boley et 
al. (2007c), such coupling can affect disk evolution. In addition, the emphasis on 
the vertical direction is well-justified for these vertically thin systems, and capturing 
the vertical transport well should guarantee that the algorithm calculates reasonable 
cooling rates. 

Consider now some incoming intensity /_ and some outgoing intensity /+. In the 
context of the approximation outlined above, the vertical flux at any cell face can be 
evaluated by computing the outgoing and incoming rays for a given column and by 
relating them to the flux with 



Once the vertical fluxes at cell faces are known, the vertical component of the diver- 
gence of the flux can be computed for the cell center by differencing fluxes at cell 
faces. 

The outgoing ray by is computed 



where At = td — tu, td is the optical depth at the base of the cell measured along the 
ray, tu is the optical depth at the top of the cell, and I+{td) is the upward intensity 



F = 2tth (7+ - 7_) . 



(2.20) 




(2.21) 
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at the base of the cell. Because I have assumed that each column in the disk is part 
of a plane-parallel atmosphere, the optical depth along the ray can be computed by 
t = r/n. 

Similar to /+, the incoming ray solution across one cell is defined as 

7_ = /_(i„)exp(-Ai) + / S{t') ex.pit^ - t')dt' , (2.22) 

where I-{tu) is the incoming intensity at the top of the cell. 

The 0th approximation for S{t) is that it is constant over the entire cell. This 
approximation leads to 

7+ = U{td) exp (-Ai) + S'o(l - exp(-Ai)) (2.23) 
7_ = 7_(i„) exp (-At) + ^0(1 - exp(-At)), (2.24) 

and 5*0 = ctTq/tt, where Tq is the temperature at the cell center. 

Because the source function is in principle a function of optical depth, additional 
complexity is necessary to obtain good accuracy. Consider a source function that 
may be represented by the quadratic 

S{t) = c + ht + at^. (2.25) 



To find the constants c, 6, and a, Taylor expand the source function about the optical 
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depth defined at the cell center to: 



S{t) ^ \ So 



dS_ 
~dt 



to 



2dt^ 



to 



dS_ 
~dt 



to 



dt^ 



to 



to >t + 



\2dt^ 



to 



e. (2.26) 



The first term in curly brackets in equation (2.26) is c, the second is 6, and the third 
is a. Using equation (2.26), I can find solutions for equations (2.21) and (2.22) across 
any given cell (see also Heinemann et al. 2006). However, in order to use equation 
(2.26), the source function's derivatives must be evaluated: 



dS 
~dt 



S' 



{T-i - T+i) 



to 

dS' 

~dt 



to 



npoKoAz 
2poKoA.z 



, and 



(2.27) 
(2.28) 



Here, the denotes the center of the cell of interest, the -1 denotes the cell center below 
the cell of interest, and the +1 denotes the cell center above the cell of interest. This 
difference scheme is used unless the following conditions are met: (A) If the +1 cell's 
density is below the cutoff value, i.e., the minimum density at which radiative physics 
is still computed, or the -1 cell's density is below the cutoff value, the derivatives are 
set to zero, which reduces the solutions for /+ and /_ to equations (2.23) and (2.24). 
(B) If cell is the midplane cell, i.e., the first cell in the upper plane, a five-point 
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~dt 



to 



,3 8To - 7T+1 - T+2 



(2.29) 



unless exception (A) is met. The simple form of equation (2.29) is due to the reflection 
symmetry about the midplane that is built into the grid, which means that the -1 
cell's values are equal to the midplane cell's values and that the -2 cell's values are 
equal to the +1 cell's values. In addition, the second derivative of the source function 
at the midplane is taken to be the average of the three-point centered difference 
method and a forward difference method, i.e., equation (2.28) is used as one would 
normally use it to compute the second derivative but that answer is averaged with 
the derivative obtained by differencing S'q and S'_^_i. Various differencing schemes have 
been tested, and this differencing scheme yields the best results for the widest range 
of optical depths and cell resolution. (C) If \AtdS' /dt\o > \2dS/dt\o, then the second 
derivative is set to zero. Condition C may appear to be strange, because it throws 
away the second derivative as soon as it becomes as important as the first derivative. 
However, this contradicts the assumption that a second-order expansion can describe 
the source function. When the second derivative is important, including only the 
first and second derivatives is inaccurate and the routine may become numerically 
unstable. I find that including the second derivative in highly optically thick disks 
along strong shocks can lead to unphysically large radiative heating, which results in 
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rapid gas expansion and disk destruction (see Figure 2.5 in §2.4.3). Condition C is not 
used for the simulation presented in Boley et al. (2007c) because the temperatures in 
the disk and the midplane optical depths were low enough that the second derivative 
always served as a correction term. 

Now that a solution for the source function integral is known, the incoming and 
outgoing intensities can be computed. The incoming ray is computed first by summing 
the solutions to the source function integral as one moves down into the disk along 
the ray with the previous sum serving as /-(t«). If desired, an incident intensity at 
t = 0, as in Cai et al. (2006), can be added to the solution by extincting the intensity 
according to the optical depth. Because reflection symmetry is assumed about the 
midplane, the incoming intensity solution at the midplane serves as the I+{tci) for the 
outgoing intensity at the midplane. 

For the r and directions, the flux-limited diffusion scheme described by Mejia 
(2004) and Boley et al. (2007c) is employed when the following conditions are met: 
(A) The vertical Rosseland mean optical depth at the center of the cell of interest is 
greater than or equal to l/-\/3. This condition ensures that we only compute flux- 
limited diffusion where photons moving vertically have less than about a 50% chance 
of escaping. (B) The cells neighboring the cell of interest also have a r > l/\/3. This 
should ensure that the code only calculates temperature gradients between relevant 
cells; the flux at this cell face is accounted for in the total energy loss (gain) of the 
system. If a neighboring cell has a r < l/-\/3, then the flux at that face is taken 
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to be the vertical flux through the flrst cell that is below r < l/-\/3 in the column 
of interest. These conditions are similar to those employed by Mejia (2004), Cai et 
al. (2006), and Boley et al. (2007c). Once fluxes have been determined for all cell 
faces, the divergence of the flux can be calculated with 



V.F^^^^^ + f^. (2.30) 

ror ro(f) oz 



2.4.3 Limiters 

There is a considerable difficulty with combining hydrodynamics, shock heating, 
and radiative heating and cooling together in an explicit scheme, namely the disparate 
timescales in concert with coarse resolution and numerical derivatives. Consider the 
hydrodynamics timescale Ath, the shock heating timescale Atg, and the radiation 
timescale Atr- The hydrodynamics timescale is the well-known Courant condition: 
th < mingridfAa;,/ (vi + Cg + cf^)}, where Axi is the cell size in direction i, Vi is ith 
velocity component of the gas, Cg is the sound speed, and cf^ — 4 (Qu)^^'^ is the AV 
diffusion speed (Pickett 1995). I define the shock heating and radiation timescales as 
Atg — fe/V and Atr = /e/ | A |, where / is some number < 1, F is the shock heating 
rate, and A is the radiation cooling/heating rate. Even though the AV timescale is 
effectively accounted for in the definition of th, strictly adhering to this definition can 
result in time steps that become extremely small and computationally inhibitive. As 
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described below, a heating limiter is employed to avoid this behavior. Likewise, a 
time step determined by tr can become too small to explicitly evolve the simulation, 
and so a radiation cooling/heating limiter is enforced. In the SV, the limiters are set 
such that lAmaxI = iTmaxI = /^/l orp. For the simulations presented by Mejia (2004) 
and Boley et al. (2006, 2007c), / = 0.1, and for the simulations presented by Cai et 
al. (2006, 2007) and Cai (2006), / = 0.03. Boley et al. (2007c) monitored the number 
of cells affected by these limiters during the calculation, and found that during the 
asymptotic phase, less than a few percent of the relevant AV heated cells were limited 
and less than a percent of the relevant radiatively cooling cells were limited. 

Although limiting the heating and cooling according to an absolute timescale 
worked for simulations presented by Boley et al. (2006) and Cai et al. (2006), which 
used the same initial model, my highly optically thick disk simulations (see Chapters 
4 and 7) often resulted in a numerical runaway. Consider the following situation. 
Over one step, too much energy, relative to an appropriate Atg, may be deposited 
into a shock and create an unrealistic temperature gradient. The following hydrody- 
namic step may be much longer than the corresponding radiation timescale, and the 
very hot cell will deposit large amounts of energy into the surrounding medium. If 
the timescales are largely disparate, this could quickly lead to a numerical runaway 
and result in the rapid expansion of gas and the destruction of the disk if an appro- 
priate limiter is not chosen (see Figure 2.5). I found that this numerical catastrophe 
was largely due to the absolute limiting timescale employed in previous simulations. 
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which was typically too aggressive in some parts of the simulation and not aggressive 
enough in others. As a result, I employ limiters that are more commensurate with my 
definitions for Atg and Atr described above. In CHYMERA, the limiters are now set 
so that heating, whether shock or radiation, can only change the internal energy by a 
small percentage, typically a percent, for every iteration on e. However, the radiative 
cooling limiter is adjusted to allow a cell to lose half of its energy for any one itera- 
tion. This limiting scheme usually leads to better stability, and it is biased toward 
faster cooling. Unfortunately, this too does not always stabilize a simulation against 
numerical runaways. I propose several possible code improvements in Chapter 8 to 
address this issue. 
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Figure 2.5 Result of a numerical radiation runaway, where radiative heating resulted 
in an unphysical hot bubble that rapidly expanded and destroyed the disk. Such a 
runaway typically occurs near strong shocks, and is believed to be caused by highly 
disparate hydrodynamics, shock heating, and radiative transfer timescales. The axes 
are in AU and the logarithmic color scale shows the midplane density and a meridional 
density slice in code units. 
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2.4.4 Opacities 

The opacities that are used in the SV and in CHYMERA are the D'Alessio et 
al. (2001) opacities. The details of the opacity tables are also discussed in Mejfa (2004) 
and in Boley et al. (2006; see Appendix A). In addition to the opacity tables, the 
SV uses the corresponding D'Alessio et al. (2001) mean molecular weight (//) tables. 
Due to an error with the inclusion of He, the typical /j, in the simulation presented by 
Mejia (2004) and Boley et al. (2006) is 2.7 instead of the standard = 2.3 for solar 
metalhcity. Inasmuch as the simulations of Cai et al. (2006, 2007), Cai (2006), and 
Boley et al. (2007c) were, in part, companion studies for the simulation presented by 
Mejia (2004) and Boley et al. (2006), the incorrect n tables were used. CHYMERA, 
in contrast, does not use the D'Alessio tables because the hydrodynamics is tied 
to the choice of In this respect, the hydrodynamics and the opacities are slightly 
inconsistent, but because the typical /i — 2.3 for the opacities, the problem is very 
minor. 

In the study presented by Boley et al. (2007c), which employed the BDNL al- 
gorithm as part of the SV, only the Rosseland optical depths were used, with the 
opacity evaluated at the cell's local temperature. This was a step backwards from 
the M2004 scheme, which employs Planck means for regions where the Rosseland 
T < 2/3. Simply switching opacities for different regions of the disk can lead to erro- 
neous physics when tracing rays in the BDNL scheme, e.g., changing the location of 
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the photosphere of the disk. Regardless, as demonstrated in §4 of Boley et al. (2007c), 
the BDNL scheme performs better overall. For CHYMERA, the mean opacity issue is 
obviated by using a midplane-weighted average between Rosseland and Planck mean 
opacities. Consider the vertically integrated midplane Rosseland optical depth r^. 
For any column of the disk, the weighted optical depth, be found by 

r.W = ^■'^ + (2.31) 

where kr and Kp are the Rosseland and Planck mean opacities, respectively. This 
method smoothly interpolates between the Rosseland and Planck regimes and ensures 
that the photosphere remains at tr = l/\/3 when the midplane optical depth is large. 
However, other techniques for smoothly switching between the Planck mean and the 
Rosseland mean opacities during the optical integration are being explored. 



2.5 Analysis Tools 

2.5.1 Fluid Element Tracer 

The hydrodynamics schemes in the SV and CHYMERA are Eulerian, and do 
not give direct information on the histories of fluid elements. In order to derive 
detailed and statistical shock information and to capture the complex gas motions in 
unstable disk simulations, I have combined a tri-Akima spline interpolation algorithm 
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with a fourth-order Runge-Kutta integrator (e.g., Press 1986) for tracing a large 
sample of fluid elements. During the integration, the Runge-Kutta scheme calls the 
interpolator each time updated velocities are required, and thermodynamic quantities, 
e.g., temperature and density, are found through interpolation at the beginning of each 
time step. Because the hydrodynamics code explicitly solves the equation of motion 
for the gas, the algorithm only needs time and velocity information to advance the 
fluid elements. 

An Akima spline is similar to a natural cubic spline, but typically yields better 
results for curves with sudden changes (see Akima 1970), as are expected in shock 
profiles. Although the spline fits a curve to a one- dimensional set of data points, 
the interpolation can be extended to data in three dimensions. First, consider a 
cubic volume with data at the vertices. A value anywhere in the volume can be 
approximated by seven linear interpolations: four to calculate values between each 
vertex in a particular direction, two to calculate the values along the projections of 
the desired point onto the interpolated lines, and a final interpolation through the 
point of interest. This scheme is depicted graphically in Figure 2.6. Extending this 
from a simple tri-linear fit to a tri- Akima spline is relatively straightforward. Instead 
of using the eight nearest points that enclose a volume, one uses the 125 closest points, 
with five data points used for every Akima spline fit. The central data point is the 
closest data value to the point of interest. I use the GNU Scientific Library Akima 
spline algorithm for performing fits. 
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Figure 2.6 Graphical representation of a tri-linear fit. The smiley face indicates the 
point for which the value is desired. The black dots indicate the vertices, the points 
for which the values are known. The values along lines adjoining the vertices are 
calculated first. The values for the projection of the smiley face onto the dark lines 
(gray circles) are then used to find the values along the hght line. Finally, the values 
for the projection of the smiley face onto the light line (white circles) are used to 
interpolate a value at the position of the smiley face. The same principle is used for 
a tri-Akima spline interpolation, but with 125 points instead of eight. 
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The integration of the fluid elements can be done by two methods: a post-analysis 
evolution or a real time integration, i.e., during the hydrodynamics evolution. For 
the post-analysis version, restart files from either CHYMERA or the SV are read 
in, typically with 1/100-1/20 of an outer orbit resolution, and the velocity field and 
thermodynamic properties are linearly interpolated between the data snapshots. This 
method is typically dominated by read-in time, and interpolation between data snap- 
shots results in smeared thermodynamic quantites and ragged profiles (e.g., density 
vs. time). The data storage demands for having sufficient data for interpolation 
are quite restrictive. However, the scheme's significant advantage is that it is post- 
analysis; different fluid element sampling with varying numbers of fluid elements can 
be used. The real time integration results in very smooth profiles, captures shocks 
well, and is not subject to inhibitive storage demands. The disadvantage is that the 
entire simulation must be rerun if a different fluid element sampling is desired. 
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Figure 2.7 Top panel: The trajectory of a fluid element in a rapidly varying velocity 
field. For comparison, the expected trajectory is also shown. Bottom panel: The 
fractional difference between the integrated trajectory and the expected trajectory is 
shown. 
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The accuracy of the interpolator and the integrator has been assessed through a 
series of tests, and is found to be satisfactory for the purpose of the fluid element 
tracer, which is to capture shock strengths and thermodynamic fluctuations. For 
example, one test is to place a fluid element in a prescribed velocity fleld: Vr — 0, 
= (GMstar/T^)^^^, and Vz = 2; sin (120). The vertical motion of a fluid element is 
then z{t) — 2;oexp (— cos (12 [Qi + 0o]) /12Q). Figure 2.7 shows the calculated and 
expected trajectories. The agreement between the known and interpolated values 
is generally within two percent. Other interpolation methods were tested, including 
a tri-cubic interpolation, a tri-cubic spline interpolation, and an inverse distance" 
interpolation, where n = 1 and 2 were tested. Generally, none of these performs 
better than the Akima spline method, although differences between the Akima spline, 
the cubic spline, and the cubic methods are slight for smoothly varying functions. In 
addition, all methods seem to be able to capture the motion of a fluid element in a 
simple velocity fleld well. Because I am interested in capturing a rapid change in the 
velocity fleld, I choose to use the Akima spline. 

2.5.2 Torque and Effective a Analysis 

The integrated gravitational torque on the inner disk due to the outer disk can be 
calculated directly from the mass distribution by 




(2.32) 
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for gravitational potential position vector x, and volume V inside r — R. When 
the integration is carried over the entire volume containing the disk, C = if angular 
momentum is conserved. As shown by Lynden-Bell & Kalnajs (1972), a symmetric 
stress tensor can be defined, where — {4:7rG)~^gigj — g'^I{87rG)~^ and dT^j/Oxj — 
—pdi^ = —pQi- Note that I use Cartesian coordinates to avoid the need to make 
distinctions between covariant and contravariant tensors. The stress tensor is related 
to the total torque on the inner disk due to the outer disk by 



where S is some surface containing the inner disk. It can be shown, as follows, that 
equations (2.32) and (2.33) are equivalent. The divergence theorem states that for 
some rank 2 tensor. 




(2.33) 




(2.34) 



Because e. 



ijkXjTki is a rank 2 tensor. 




(2.35) 



Recall that dxj/dxi = 5ji by definition, and so 




(2.36) 
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Inasmuch as Tj^ is a symmetric tensor, eij^Thj — 0, and by definition of 7}^, 



which is equation (2.32). 

In the analysis of these three-dimensional disk simulations, I use equation (2.32) 
to calculate the gravitational torque because the volume intergral can be neglected in 
regions where the density goes to zero. Moreover, I am only interested in the torque 
in the direction, and so the torque intergral can be rewritten as 



Using this description for the torque at any given radius, one expects for extrema in 
the torque profile to indicate transitions between mass inflow and outflow. 

In addition to the gravitational torque, velocity fluctuations in a mean flow can 
create Reynolds stresses (Landau & Lifshitz 1987). The torque in the tz direction 
due to these fluctuations for a given r is given by 



where the Reynolds stress tensor for the r-<p component is Tj.^ = p6vrSv^. The 
fluctuating components of the velocity fleld Svi are found by subtracting the mean 




(2.37) 




(2.38) 




(2.39) 
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flow from the velocity at a given cell. Calculating an accurate Reynolds stress is 
therefore dependent on determining a meaningful average flow for the gas. In the 
analysis presented here, I flnd the mean flow by calculating density weighted velocity 
averages for each r in the disk. The Reynolds stresses determined this way are still 
extremely noisy, and it is unclear whether they represent the actual Reynolds stress 
contribution to the total torque well. This difficulty is due to the relatively coarse 
resolution used in these global simulations. As will be shown in Chapter 6, the 
gravitational torque alone represents the mass transport in these disk simulations 
fairly well, and so I typically ignore the Reynolds stress. 

Once the torque on the disk is determined, an effective a can be calculated, i.e., 
the Shakura & Sunyaev a parameter needed to describe mass transport if the disk 
were an a disk, where (Gammie 2001) 



Here, % is the vertically integrated, azimuthally averaged total (gravitational+Reynolds) 
stress tensor and Cg is the sound speed deflned for a razor-thin disk. Gammie (2001) 
showed that if the disk evolved toward a state where local disk cooling was balanced 
by local shock heating and gravitational work, the effective a can be described by 



(2.40) 



a = 



d In r c^E 




(2.41) 
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for a Keplerian disk, where ^cooi is the local cooling time and 72D is the two-dimensional 
adiabatic index. The two- and three-dimensional 7s are related by 72D — 3 — 2/7 and 
72D — (37 — 1) / (7+ 1) for the strong and negligible self-gravitating limits, respectively. 
Because the torque in the vertical direction is desired, the surface on which to evaluate 
% is an infinite cylinder centered along the z axis. If only the gravitational stress 
tensor needs to be considered, then the azimuthally averaged, vertically integrated 
stress can be related to equation (2.32) by 



In order to calculate effective a profiles for the simulations presented in this disser- 
tation, I use midplane sound speeds and equation (2.41) to evaluate equation (2.40). 
Because the vertically integrated in equation (2.40) is probably less than the mid- 
plane c^, this approach likely understimates a. 




(2.42) 



Chapter 3 



RADIATION TESTS 



A crucial demand on any radiation algorithm designed for the protoplanetary disk 
problem is the necessity to handle both the high and low optical depth regimes. An 
analytic solution to a relevant test problem must be found in order to evaluate the 
accuracy of radiative transport algorithms in disks. In this Chapter, I describe a toy 
problem based on Hubeny (1990) that can be used to test the accuracy of a radiation 
routine for a disk geometry. 

Consider a plane-parallel slab with constant vertical gravitational acceleration g 
but with a midplane about which reflection symmetry is assumed. Suppose there 
is some heating mechanism that produces a known distribution of astrophysical flux 
once the system reaches hydrostatic and radiative equilibrium; in equilibrium, energy 
transport is only vertical. Make the ansatz that the vertical astrophysical flux has 
the form 



where To — ctT^/t^ is the astrophysical flux from the atmosphere with effective tem- 




(3.1) 
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perature Tg, r is the Rosseland mean optical depth measured vertically downward, 
and Tm is the optical depth at the midplane. This function ensures that the flux goes 
to zero at the midplane and that J^^ — t — 0- The heating term required to 

achieve this flux distribution is then 



r = -.^ = .^A (3.2) 



where p is the density at the point of interest and k is the Rosseland mean mass 
absorption coefficient. 

If '^m ^ 10, the temperature structure may be derived from the flux by using 
the standard Eddington approximation, which relates the mean intensity J to the 
astrophysical flux by 

4 dJ (r) 

3^ = ^.(r). (3.3) 

Integrating equation (3.3) yields 



1-^ 
2t^ 



+ q). (3.4) 



The constant q can be determined by considering the low optical depth limit. In that 
limit, the atmosphere reduces to a standard stellar atmosphere; therefore q = l/\/3 
(Mihalas & Weibel-Mihalas 1986). 

In the limit that Tm ^ 0.5, the atmosphere approaches an isothermal structure. 
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Because the source function becomes constant, the observed flux can be found from 
F — 27r/x/+ where /j, = l/\/3 for the single ray scheme (Chapter 2), which gives the 
temperature of the isothermal atmosphere: 

2/i(l-exp[-2W/i])' ^""-^^ 

where the factor of 2 in the exponential accounts for both sides of the atmosphere. In 
addition, 1 explicitly include because the vertically integrated r is used. As <^ 1, 

An additional assumption about the form of the opacity law permits analytic eval- 
uation of the hydrostatic structure of the disk and allows one to control whether the 
atmosphere will be convective. For example, I assume that k is constant throughout 
the disk for two of the tests presented below. This assumption should make the disk 
convectively stable (Lin & Papaloizou 1980; Ruden & Pollack 1991). It also makes 
the relation between pressure and optical depth simple: 



p = —T. (3.6) 

K 



Moreover, the midplane pressure is related to the full disk surface density S by 



Pm = 9^- 



(3.7) 
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3.1 Relaxation Test 

For the relaxation test, I heat the gas slab in accordance with equation (3.2) 
to test whether the radiation routine in question allows it to relax to the correct 
hydrostatic and radiative equilibrium configuration. By selecting different midplane 
optical depths, the effects of resolution on the routine can be tested. 

Following Boley et al. (2007c), I test the radiation routines used by M2004, C2006, 
and BDNL. The initial condition for the atmosphere is a constant density structure 
with a, — 0.05, 0.5, 5, 10, and 100. The target effective temperature of the 
atmosphere, which controls the magnitude of the heating term, is set to Tg = 100 K, 
and k — Kq— constant. For the discussion that follows, I take the high optical depth 
regime to indicate where r > 1/a/3 and the low optical depth regime to indicate 
where r < l/-\/3. Note that this boundary is slightly different from the M2004 and 
C2006 schemes; these schemes use r = 2/3 to set the boundary between the low and 
high optical depth regimes, as is done in the standard Eddington solution. 

Figure 3.1 compares the M2004, C2006, and BDNL solutions to the analytic tem- 
perature profile and to the fiux profile for = 100. The BDNL routine matches the 
analytic curves very well. The M2004 routine does well in matching the temperature 
curve for most of the high optical depth regime, which is resolved by 12 cells, but 
the low optical depth regime has a sudden temperature drop. As reported by Boley 
et al. (2006), this temperature drop is an artifact of the lack of complete cell-to-cell 
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coupling in the optically thin region. Despite the temperature drop, the boundary 
condition between the two regimes typically yields the correct boundary flux, i.e., the 
flux leaving the high optical depth regime. An unfortunate result of this boundary 
condition is that it produces oscillations in the flux profile and in the temperature 
profile when the low/high optical depth boundary is near a cell boundary or when the 
entire disk height is only resolved by a few vertical cells. These numerical oscillations 
produce artificial heating, and it is this behavior that may contribute to keeping the 
inner 7 AU of the disk presented in Boley et al. (2006) hot. Although this is problem- 
atic, even when oscillations occur, the time-averaged cooling time is within 10% of 
the expected value for the test shown in Figure 3.1. The C2006 improvement lessens 
this problem. 

Figure 3.2 compares the results of the M2004 and BDNL routines with each other 
for Tm — 10. The high optical depth regime is now resolved by only six cells. Both 
methods compute the correct flux through the slab. The temperature profiles are both 
skewed more than in the Tm = 100 case mainly because the solution deviates from the 
Eddington approximation, which is used to derive the analytic temperature profile, 
as Tm becomes small. Figure 3.3 shows the same comparison as done in Figures 3.1 
and 3.2, but for = 5. The high optical depth regime is now resolved by only four 
cells. Again, both methods allow for the correct flux through the slab, and the BDNL 
temperature distribution is close to the analytic value, with slight departures again 
due mainly to the inaccuracy of the analytic curve. 
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To demonstrate the accuracy of each algorithm in the low optical depth limit, 
I show the temperature profiles in Figure 3.4 for the M2004 routine and the BDNL 
routine for = 0.05 and 0.5, and compare each curve with the temperature estimate 
calculated from equation (3.5). The M2004 routine yields temperature profiles that 
are colder than the expected temperature profiles, and the departure is more severe 
as increases. This is a result of the lack of complete cell-to-cell coupling in the 
atmosphere. The BDNL routine is in excellent agreement with the predicted temper- 
ature. The departure from the analytic estimate observed in the = 0.5 case is a 
result of the inaccurate assumption in the analytic curve that the source function is 
constant for r ~ 1. 
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Figure 3.1 Results of the relaxation test for the M2004, C2006, and BDNL schemes 
with Tm = 100. Top: The left panel shows the relaxed temperature profile for the 
M2004 and C2006 routines, while the right panel shows the relaxed temperature 
profile for the BDNL routine. In both panels, the analytic curve is represented by the 
red, dashed curve. The first sudden drop in the M2004 and C2006 schemes is due to 
the lack of complete cell-to-cell coupling required by radiative transfer. The second 
drop, which is also in the BDNL profile, is where the density drops to background and 
where the algorithm stops following radiation physics. Bottom: Similar to the top 
panels but for the flux profile. The undulations in M2004's flux profile are believed 
to be due to the low/high optical depth boundary lying near a cell face intersection. 
The C2006 modification helps avoid this problem. Finally, the sudden drop in the 
M2004 and C2006 flux profile occurs because these schemes only explicitly track the 
flux through the optically thick disk. The right panels are the same as Figure 1 in 
Cai et al. (2007). 
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Figure 3.2 Same as Figure 3.1, but with = 10 and only for the M2004 and BDNL 
schemes. The undulations in the M2004 flux profile are no longer present. The slight 
departure of the BDNL solution from the analytic temperature curve is mainly due 
to the breakdown in the Eddington approximation, which is assumed in the analytic 
temperature profile, as becomes small. 
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Figure 3.3 The same as Figure 3.2, but with = 5. Once again, the departure of the 
BDNL temperature solution is mainly due to the breakdown of the approximations 
used to calculate the analytic curve. 
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Figure 3.4 Temperature curves for the low optical depth limit. The red, dashed curve 
indicates estimated isothermal temperatures from equation (3.5), the solid, dark curve 
indicates the results of the BDNL routine, and the dashed, dark curve indicates the 
results of the M2004 routine. The curves at the lower temperature correspond to 
Tm = 0.5, and the curves at the higher temperature correspond to Tm = 0.05. The 
lower optical depth corresponds to the higher temperature because cooling is less 
efficient. The M2004 routine yields temperatures that are too cold, because it always 
uses the free-streaming approximation when r < 2/3. However, the M2004 routine 
converges to the correct solution as — > 0. The BDNL routine is consistent with 
the estimated temperature for the Tm = 0.05 case, and it is roughly consistent with 
the Tm = 0.5 case. The inconsistency seen in the Tm = 0.5 case is a result of the small 
Tm assumption, which is used to derive equation (3.5). 
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3.2 Contraction Test 



In order to study how accurately the radiation algorithms work with the hydro- 
dynamics routines and to study the effects of resolution, I allow the slab to cool 
and follow the contraction. If one assumes a constant opacity law and a large r^, 
then the contraction becomes homologous, and a relationship between the midplane 
temperature and the cooling time is easily attainable. Consider the cooling time 



where /i ~ E/p^ is the scale height of the atmosphere and U is the internal energy 
per unit area. If one assumes an ideal gas law, then one expects that 



U 



Pmhr, 



m 




1 



(3.9) 



^cool 




Finally, because U ^ Pmh ^ Tj 



(3.10) 
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Figure 3.5 Contracting slabs as shown in the tcooi-f^~^ plane for the same cases as 
shown in Figure 3.1; all curves should be linear if the slabs contract as expected. Both 
schemes break down when the high optical depth regime is contained within 5 cells 
{U~^ pa 0.75), but the M2004 routine (light curve) starts to deviate from the expected 
solution once the slab is resolved by 6 cells {U~^ ^ 0.43). The sudden decreases in 
cooling times for the M2004 routine occur when the optically thin/thick boundary 
transitions into another cell. The dotted line shows the expected behavior. 
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For this test, I take the relaxed slab shown in Figure 3.1 with = 100 and 
turn the heating term off. The contraction is followed until the scheme breaks down. 
Figure 3.5 demonstrates this contraction sequence in the tcooi-U~^ plane when the slab 
is evolved with the BNDL scheme (heavy, dark curve) and when it is evolved with 
the M2004 scheme (lighter, gray curve). Both cases follow the expected contraction 
closely, which is indicated by the dotted line, until the optically thick atmosphere is 
resolved by five (BDNL) or six (M2004) cells. The sudden decreases in cooling times 
for the M2004 routine occur when the optically thin/thick boundary transitions into 
another cell. 

3.3 Convection Test 

Another test I describe demonstrates whether the radiation scheme permits con- 
vection when it should occur. Lin & Papaloizou (1980) and Ruden & Pollack (1991) 
show that convection is expected in a disk-like atmosphere when the Rosseland mean 
optical depths are large and when (3 > /3crit for k ^ T^. For a 7 = 5/3 gas, the 
vertical entropy gradient is driven to a negative value when the critical exponent 
/Scrit ~ 1-5. Thus, I present two cases for which almost identical atmospheres are 
allowed to relax to an equilibrium. In one case, P — 1, which should make the atmo- 
sphere convectively stable, and the other has P — 2, which should make it unstable. 
As shown in the top panels of Figure 3.6, I find that the BDNL routine produces 
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convection when it should and does not when it should not. Likewise, the bottom 
panels of Figure 3.6 demonstrate that the M2004 routine also permits or inhibits 
convection correctly. However, M2004 does seed an artificial superadiabatic gradient 
at the boundary between the optically thin and thick regions due to the temperature 
drop at the photosphere. Nevertheless, convection does not occur for /3 = 1 even 
with this seed, and the superadiabatic gradients are an order of magnitude smaller 
for /3 = 1 than for /3 = 2 at that boundary. 

Finally, I measure the flux that is carried by convection in the (3 — 2 case for 
the BDNL scheme. As indicated in Figure 3.7, convection carries about 30% of the 
flux. Figure 3.7 also acts as a reminder that the energy must ultimately be radiated 
away. Boley et al. (2006) and Raflkov (2007) argue that because the convective flux is 
controlled by the radiative flux leaving the photosphere of the disk, convection should 
not be expected to lead to rapid cooling and fragmentation in protoplanetary disks, 
as claimed by Boss (2004a) and Mayer et al. (2007) (see Chapter 5). 
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Figure 3.6 Convection test. The heavy black contour indicates the same density 
contour for each panel. Arrows show relative velocities in the r and z directions. The 
blue contours indicate superadiabatic regions. The motions in the left panels are a 
few orders of magnitude smaller than the motions in the right panels. Top: BDNL. 
The left panel shows the case /3 = 1; convection and superadiabatic gradients are 
absent, and the velocities represent low level noise. The right panel shows the case 
f3 = 2; convective cells and superadiabatic gradients are present. Bottom: Same as 
the top panel but for the M2004 scheme. The superadiabatic regions near the top 
density contour are due to the artificial, sudden drop in temperature at the optically 
thin/thick interface. The superadiabatic gradients in the left panel are about an order 
of magnitude smaller than those in the right panel. 
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Figure 3.7 Average flux through the atmosphere with (3 = 2 and with the BDNL 
scheme. Convection carries about 30% of the flux at maximum, but almost all the 
flux is radiative in the photospheric region near cell 15. 
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3.4 2D Test 

The M2004, C2006, and BDNL radiation routines include more cell-to-cell cou- 
pling in the vertical direction than the radial or azimuthal directions. Vertically 
integrated optical depths are also used to determine the radiation boundary condi- 
tions. Generally, this vertical bias is justified in the highly fiattened disk systems that 
are evolved in the SV and in CHYMERA. However, the degree to which temperature 
structures are skewed by this bias should be investigated using 2D and 3D tests. Al- 
though the 2D test is still being developed, I present a few results here for the BDNL 
routine. 

A spherical, constant density mass distribution is loaded onto the grid and held 
fixed. This isolates the radiation algorithm in the code. A constant energy input 
source is included in just the center cells, i.e., the first computational cell closest to 
the r, z origin for all (f). The opacity is set to some constant, and the magnitude of 
the constant is adjusted to vary the optical depth regimes for different tests. Figure 
3.8 shows the results. The temperature contours are roughly spherical for the high 
optical depth case, which indicates that the optically thick/thin boundary conditions 
are reasonable. 

The low optical depth case does not perform nearly as well. For this regime, 
the temperature contours are highly skewed along the r axis, and there is an odd 
temperature inversion at the photosphere near the z axis (Fig. 3.8). At optical depths 
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of about 5, the radiation transport is not quite in the diffusion hmit (see §3.1), and 
so fluxes along and are disparate. This creates the temperature skew along the 
r axis because radiation from the heat source in the center cell is coupled over several 
cells in the z direction by the ray solution, but is only coupled between adjacent cells 
in the r direction by the diffusion equation. 

Near the z axis, the temperature inversion in the photosphere is also due to 
the difference in handling radial and vertical fluxes. The combination between the 
flux-limited diffusion approximation in low-optical depth regimes and the radiation 
boundary condition described in §2.4.2 introduces too much flux along r, and these 
boundary cells cool more than expected. When the flux along the r direction is set 
to zero for the boundary cells, the temperature inversion is corrected. 

In addition to temperature contours, the fluxes along the r and z axes are shown 
in Figure 3.9 for the high and low optical depth cases. The fluxes for the high optical 
depth case follow the expected x"'^ profile. There is some deviation near r — z — 0, 
but this is likely due to the cylindrical energy input source near the center of the 
sphere. The low optical depth case shows much more deviation from the expected 
profile. As discussed in Chapter 6, routines that include complete cell-to-cell coupling 
in all directions is a critical next step for radiation algorithms in disk codes. Although 
it should be noted that for both optical depths, the total luminosity of the sphere 
reaches the expected value. 
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Figure 3.8 Temperature (left panels) and vertically integrated r (right panels) plots 
for the 2D test of the BDNL algorithm. The high optical depth test (top) shows 
a slight flattening of the temperature contours, but is satisfactory overall. The low 
optical depth test (bottom) shows much more flattening of the temperature contours 
and odd behavior along the r axis and near the z axis at the photosphere (see text). 
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Figure 3.9 Vertical and radial flux profiles for the 2D high (top) and low optical 
depth (bottom) tests. The high optical depth profiles converge to an x~'^ profile, as 
expected. There is a deviation (only shown in part) as one approaches r — z — 0, 
but this may be due to the nonspherical energy source. The low optical depth case 
shows odd behavior. For the flux along the z axis, the flux strongly deviates (only 
shown in part) from the expected profile near r = z = 0, but converges as one moves 
outward. In contrast, the fiux along the r axis follows an profile near the center, 
but deviates as one moves outward. The kinks at the very end of each fiux profile 
are due to the edge of the sphere. Although not understood at the moment, the 
magnitudes of the fiux profiles are also different. 
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INITIAL MODEL GENERATION 

In this Chapter I describe the various initial models that are used in these studies. 

4.1 The Mejia, Light, and Moderate Disks 

As discussed in Mejia (2004), Mejia et al. (2005), and Cai (2006), initial conditions 
can be generated as star+disk equilibrium models using a modified Haschisu (1986) 
self-consitent field (SCF) relaxation method. In this method, the density and angu- 
lar momentum distributions, assuming a poly tropic equation of state, are iteratively 
solved until convergence to the desired model is reached. A star+disk model is created 
by specifying Mstar/^totah ^poi/^O) and p, where Tq is the outer, initial edge of the disk, 
Tpoi is the polar radius of the protostar, and E(r) ~ is the desired surface density. 
During the iterations, the angular momentum distribution is adjusted to create the 
desired law for the disk, while the protostar is specified to be either a uniformly ro- 
tating or non-rotating polytrope. The goodness of the equilibrium model is evaluated 

by measuring the virial test of the model: VT = |ii^grav + 2-E'rot + 3P | / 1 -Egrav | , where 
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V — J pdV and VT — for perfect equilibrium. 

Because the star+disk model is created in dimensionless code units, the model 
is scalable. However, scahng the disk beyond a few AU results in an unreahstic size 
for the central star. In addition, the hydrodynamics time step is restricted by the 
Courant condition, and so the hot star severely limits the time step size, which is 
problematic for long disk simulations. To obviate these difficulties, the central star 
is removed from the grid, but its gravitational potential is stored (see Mejfa 2004) or 
replaced by a point mass potential with the origin at the center of the grid. Replacing 
the potential with a point mass only creates small disturbances in the equilibrium 
model, and it allows for easy expansion of the grid. For implications regarding keeping 
the central mass' potential fixed, see Chapter 6. 

For the studies presented here, three of the equilibrium models were generated 
by the SCF method: the Mejfa disk, the moderate disk, and the light disk. Their 
model parameters are given in Table 3.1. The Mejia disk has been used for the initial 
conditions in multiple studies, including Mejia (2004), Mejia et al. (2005), Boley et 
al. (2005), Cai et al. (2006, 2007), Cai (2006), Boley & Durisen (2006), Boley et 
al. (2006, 2007a, 2007c). The minimum Q for the Mejfa disk is about 2. The light 
disk and moderate disks were used as the initial conditions for the shock bore studies 
(see Chapter 5) discussed in Boley et al. (2005) and Boley & Durisen (2006). Both 
disks are stable against GIs, with the minimum Q ^ 10 for the light disk and 5 for 
the moderate disk. For Boley & Durisen (2006), the moderate disk was cooled in a 
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2D version of the SV (Pickett 1995; Mejia 2004) so that the minimum Q 2. Figure 
4.1 shows density contour plots scaled to the central density of the now removed 
protostar, p^, for each model. 



The flat-Q disk is created for evolution with a variable first adiabatic index with an 
ortho/para = 3:1 (Chapter 2). To create an equilibrium model with the SCF method, 
significant modifications would be required. Moreover, even with a simple adiabatic 
index, it was very difficult for the SCF method to converge to an equilibrium solution 
for the desired model parameters (below) . Because there was little time for extensive 
code modifications and because computer cycles were readily available on department 
rendering machines, a far less eloquent method was used. 

Consider a disk that is vertically polytropic, is axisymmetric, has a constant 7, 
and is roughly Keplerian. Also assume that self-gravity is negligible for the vertical 
structure. For this case, I find that 



where K is the polytropic coefficient for any given r, G is the graviational constant. 



4.2 Flat-Q Disk 




(4.1) 
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and Q, is the Keplerian angular speed. To create a model, one needs to specify power 
laws for E and Q, E and K, or K and Q. In practice, I specify E in order to control 
the mass of the disk and then choose either a constant K profile or a constant Q 
profile. Once the profiles have been determined, the disk structure can be derived by 
using the following equations. 

p{z) = po{l-zyh'Y/^'-'\ (4.2) 
= f]2(^-l) ^o"'' ('^■3) 

- 1 [2^) r[7/(7-i)] ' ' ^'-'^ 

where h is the disk scale height. Note that when the K is constant and Q is constant, 
the surface density profile follows a r~^-^'y profile. 

Equations (4.1-4.4) allow one to solve for a model disk. However, the disk will 
not be in equilibrium because a Keplerian rotation was assumed and self-gravity was 
neglected. Moreover, a constant 7 was still assumed. Why use this analytic method 
over the SCF method? (1) Both methods are guaranteed to produce a model that 
is out of equilibrium when loaded onto the grid if a variable 7 is introduced. (2) 
The SCF method has difficulty converging models for disk parameters desired for the 
flat-Q case. (3) The model is intended for a radiation hydrodynamics simulation, 
and although both the SCF method and the analytic model produce a polytropic 
model in rough equilibrium, the models are out of radiative equilbrium. Because 
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both methods require some relaxation in the hydrodynamics code, the easiest and 
most straightforward method is the analytic model. Unfortunately, the relaxation 
of the analytic model to an equilibrium structure proves to be somewhat difficult, 
and requires employing numerical tricks, e.g., alternating between dampening the 
momentum in the radial and vertical directions. Some oscillations of the outer disk 
could not be dampened, and remain in the initial model. 

By using the analytic method described above, the flat-Q model is produced with 
7 = 5/3, E ~ r~^, Q = constant, and K r. The total disk mass is approximately 
O.O9M0, and in accordance with the E and K profiles, the midplane temperature 
profile T ~ . The initial disk has a radius of 10 AU, and the surface density 
at 5 AU is 3300 g cm~^. The actual Q for the model after relaxation is about 2 
everywhere. 

In order to study spiral waves, shocks, and mass transport when a dead zone 
becomes highly unstable, the flat-Q model is modified to drive a region of the disk 
toward instability. Mass is added to the disk, linearly in time, in a Gaussian profile, 
with the peak centered on 5 AU. The internal energy density is adjusted so that 
the specific internal energy remained the same. The radial FWHM of the Gaussian 
is chosen to be 3 AU, which is roughly the size of the most unstable wavelength 
Au 27r/i 0.27rr (Durisen et al. 2007b). Mass is added until nonaxisymmetry is 
visible in midplane density images, which took approximately 190 yr, i.e., 6 orp of 
evolution (1 orp is defined as 1 outer rotation period at 10 AU for this model). The 
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total mass added is about 0.08 M©, and so if one were to imagine a corresponding 
accretion rate, it would be about 4 x 10~^ Mq yr~^. It is probably unphysical to add 
so much mass to the disk so quickly without increasing the temperature of the disk. 
However, I remind the reader that the study is a numerical experiment, as described 
in Chapter 7. 

The mass build-up drops Q below unity in a narrow region around 5 AU just before 
the disk bursts. Because the ring is grown over 6 orp, the instabilily is probably not 
overshot. The mass-weighted average Qav over the FWHM of the Gaussian centered 
at 5 AU is approximately unity when the ring becomes unstable. This implies that 
Q must approach unity for a spatial range of at least an unstable wavelength before 
GIs will activate in a disk. 

Figure 4.2 shows density contours for the flat-Q model and the flat-Q model with 
the mass concentration. As in Figure 4.1, the density is normalized to a pc, which 
is based on the typical magnitude difference between the maximum disk density and 
maximum protostar central density in the SCF models. 
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Table 4.1 Model parameters for initial disks generated by the SCF method. 



Model 


Mstar/Mtotal 


r-poi/ro 


P 


VT 


Qmin 


Mejia 


0.875 


20 


0.5 


2.6e-3 


2 


Light 


0.982 


10 


0.5 


6.5e-4 


10 


Moderate 


0.964 


10 


0.5 


6.3e-4 


5 


Cooled Moderate 


0.964 


10 


0.5 


6.3e-4 


2 
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Figure 4.1 Density contours, with p/pc — l.e-4, l.e-5, l.e-6, l.e-7, and l.e-8, for the 

Mejia disk, the hght disk, the moderate disk, and the cooled moderate disk, from the 
top left to bottom left respsectively. The odd density enhancement at the inner edge 
of the light disk is an artifact of removing the protostar; it is partly accreted to the 
protostar and partly subsumed into the disk during the beginning of the simulation. 
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Figure 4.2 Similar to Figure 4.1, but with p/pc = l.e-3, l.e-4, l.e-5, l.e-6, l.e-7, and 
l.e-8 for the flat-Q model. The left panel is the initial flat-Q model and the right 
panel is the same model, but with the density enhancement. For these models Tq = 10 
AU. 



Chapter 5 



SHOCK BORES 

5.1 Introduction 

To address the dynamic evolution of protoplanetary disks properly, the three- 
dimensionality of spiral waves and shocks must be understood. Fully three-dimensional 
waves in accretion disks have been studied by several authors (e.g., Lubow 1981; Lin 
et al. 1990; Korycansky & Pringle 1995; Lubow & Ogilvie 1998; Ogilvie 2002a,b; Bate 
2002) , typically in the context of tidal forcing. It has been noted (see Lubow & Ogilvie 
1998, hereafter L098) that these waves act like fundamental modes (/-modes), which 
correspond to large surface distortions in the disk. These waves can affect the disk's 
evolution through wave dissipation at the disk's surface and through gap formation 
(Bryden et al. 1999; Bate et al. 2003). The ability of radially propagating waves to 
transport angular momentum, and also influence gap formation, is dependent on the 
thermal stratification of the disk (e.g., Lin et al. 1990; L098). Finally, consider the 
following time scale argument. The dynamical time at a given radius in a disk is 
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about the orbital period, i.e, ta — Prot — SttQ, where is the orbital frequency. The 
vertical time scale can be thought of as the time it takes for a wave launched at the 
midplane to propagate to the disk surface and back. If h is the disk scale height, and 
c is the midplane sound speed, then the vertical time scale ty — 2h/c 2Q = Pj-qi/t:. 
To within a factor of order unity, which depends on the thermal stratification of the 
disk, the time scales are commensurate. Thus the vertical direction, although often 
ignored, plays a very important role in disk dynamics. 

Independent work in the context of gravitational instabilities (GIs) has also found 
that once GIs become well developed, e.g., in the asymptotic state of Mejia et 
al. (2005), the spiral waves resulting from GIs in the disk also behave like /-modes and 
involve large surface distortions (Pickett et al. 1998, 2000, 2003; Durisen et al. 2003). 
Furthermore, it was noted in these studies that the spiral wave activity is highly 
nonlinear and involves many modes of comparable strengths (see also Gammie 2001; 
Lodato & Rice 2004). In the asymptotic disk, sudden increases in disk scale height 
occur. These splashes of disk material seem to be related to shocks in the disk and 
could have important consequences for disk evolution through the generation of tur- 
bulence, breaking waves, and chondrule formation (Boley et al. 2005, 2006). In order 
to understand this shock-related splashing, which is evidently highly nonlinear, an 
approach other than /-mode analysis is required. 

Martos & Cox (1998, hereafter MC98) investigated shocks in the Galactic disk, 
where the otherwise isothermal equation of state (EOS) is stiffened by magnetic fields. 
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Their findings indicate that shocks occurring in semi-compressible fluid disks have 
characteristics of hydraulic jumps and behave, in part, like gravity modes (gf- modes). 
A classical hydraulic jump occurs in an incompressible fluid, where the only way 
to reduce the kinetic energy of fluid elements coming into the wave is to convert 
it to gravitational potential or turbulent energy (e.g., Massey 1978). Examples of 
hydraulic jumps can be found in spillways wherever there is an abrupt change in 
slope and the slowly moving water has a greater height than the rapidly moving 
water. The hydraulic jump is the abrupt change in the height of the flowing water. 
The same phenomenon can be exhibited in a disk because, like the flowing water, 
the unperturbed disk is usually hydrostatic in the z direction and work must be 
done against gravity to expand the disk in the vertical direction. Abrupt changes 
in the scale height of a disk are thus possible for similar reasons. MC98 present an 
analytical theory, along with two-dimensional magnetohydrodynamics simulations, 
in the context of hydraulic jump conditions, and note that the jumping gas can 
lead to high-altitude shocks when the jumping material crashes back onto the disk. 
After MC98, Gomez & Cox (2002) simulated a portion of the galactic disk in three 
dimensions, noting a behavior of the gas similar to that found by MC98. However, 
their analyses neglect self-gravity, which can affect the morphology of a global shock. 

I refer to these hybrids between hydraulic jumps and spiral shocks as shock bores. 
Following Boley & Durisen (2006), I lay out in this Chapter a theory that combines 
the Rankine-Hugoniot jump shock conditions for both the adiabatic and isothermal 
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cases with the classical hydraulic jump conditions. Fully three-dimensional hydrody- 
namics simulations are discussed in §5.3 to illustrate simple cases of spiral shocks in 
protoplanetary disks. I discuss some implications of these results in §5.4 and summa- 
rize the main conclusions of this study in §5.5. 

5.2 Shock Bores 

5.2.1 Plane- Parallel Approximation 

In a disk, abrupt vertical expansions will typically occur after a shock. The 
reason for a shock bore can be understood by considering hydrostatic equilibrium 
(HE) and the EOS. For simphcity, assume that (1) the shock is planar, (2) it is 
propagating in the x direction, (3) the disk is vertically stratified in a vertical direction 
z perpendicular to the x direction, with the pre-shock region in vertical HE, and (4), 
except for the discontinuity at the shock front, ignore variations in the x direction. 
With conditions (3) and (4) one may write 



1 dp 
p dz 



dz 



(5.1) 



for the pre-shock flow, where $ is the total gravitational potential, p is the pressure, 
and p is the density. Consider first the case of an adiabatic shock. Then for the 
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Rankine-Hugoniot shock conditions, 
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^ = ^ - 1^ and (5.2) 
Pi 7 + 1 7 + 1 



ui P2 7 + 1 7 + lA^'' ■ 

where subscripts 1 and 2 represent the pre- and post-shock regions, respectively, and 
where the gas speed relative to the shock is u, the ratio of specific heats is 7, and M. 
is the Mach number given by 

= ^. (5.4) 

7Pi 

If the pre-shock gas is initially in vertical HE, what is the state of vertical force 
balance behind the shock? For the adiabatic case, using equations (5.2) and (5.3), 
one can write the ratio of the post-shock and pre-shock pressure body forces as 



]^dp2 r ]^dpA~^ ^ / 27M^ 7-l\ 
P2 dz \pi dz J V7 + 1 7 + 1/ 



7-1^ 



7 + 1 (7 + 1) 
27M^ (7 - 1) - (1 - 67 + 7^) - 2 (7 - 1) 
M2(7 + l)' 



(5.5) 
. (5.6) 



For the gravitational body force, let be the background potential, presumably due 
to the primary, and let and be the contribution to $ from the gas self-gravity 
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in the pre- and post-shock regions. The ratio of the potential gradients then becomes 



d^2 f d^i\ dz^^ + dz^ 
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dz \ dz J c^2$* + dz^g-^ 

g + 1 



(5.7) 
(5.8) 



where 



, = ^fd?^V' (5.9) 



dz \ d^; 

is the ratio of the background and the pre-shock gas potential gradients. Condition 
(4) permits one to write ^ V^^ oc p; therefore, the ratio of the self-gravity 
potential gradients is equivalent to the ratio of the densities. This yields 



d^2 fd^i\ ^ q + a 



where 
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Define the jump-factor Jf to be the ratio of equations (5.5) and (5.7): 



Jf 



1 dp2 f 1 dpi 



P2 dz \pi dz 



d^i { d^ 



dz \ dz 



2jM^ (7 - 1) - (1 - 67 + 7^) - 2 (7 - 1) 
^2(7 + 1)' 



g + 1' 
q + a 



(5.12) 
(5.13) 



Note the hmits of equation (5.13). When self-gravity dominates {q — > 0) and A4 is 
large, 



Jf 



27M2 _ 



(7+1) 



3 ' 



(5.14) 



and, when the background potential dominates {q — > 00) and A4 is large. 



Jf 



(7 - 1) 
(7 + If 



(5.15) 



When self-gravity is negligible, Jf becomes the post- and pre-shock temperature ratio 
for a purely adiabatic shock. To understand the significance of J/, consider the vertical 
acceleration of the gas in the post-shock region: 



dvz 
dt 



1 dp2 _^ d^2 
P2 dz dz 



(5.16) 
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Using equation (5.12) in (5.16) yields 



d^2 

dz 



Jf 




+ 1 



(5.17) 



d^2 




(5.18) 



dz 



In the limit of no self-gravity, where g — > oo, ~ 



z {Jf — 1) ior a thin disk, where 



fl is the Keplerian circular angular speed. Equation (5.18) demonstrates the physical 
meaning of Jf. When Jf > 1, the gas is over pressured, and it will expand vertically, 
while for Jf < 1, self-gravity causes the gas to compress. If the shock is truly strong, 
an expansion will always occur in the adiabatic case (7 > 1). However, self-gravity 
and 7 alter the strength of the adiabatic vertical acceleration (see Fig. 5.1). For low 
7 and q, regimes with Jf < I exist for low A4. 

Now consider the case of an isothermal shock, which reduces the shock conditions 

to 



For these conditions with the same assumptions as stated for the adiabatic case, the 
resulting Jf is 




(5.19) 



Jf 



q + l 



(5.20) 
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02468 10 02468 10 



Mach Number Mach Number 

Figure 5.1 The three curves show the adiabatic Jf with respect to for g — > 
oo (soUd), 1 (dotted), and (dashed), corresponding to no self-gravity, equal back- 
ground and self-gravity, and fully self-gravitating, respectively. Notice that for 
7 = 5/3 (left), the disk only collapses for the self- gravity dominated case and only at 
low Mach numbers. However, when 7 = 7/5 (right), the compressions are noticeable 
at low Al even when self-gravity is only just comparable to the background potential. 
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What one can immediately see from equation (5.20) is that if self-gravity is important 
(typically q < 100), one expects an isothermal gas to compress vertically in the post- 
shock region because M. > 1. If ^ — > oo, spiral shocks in an isothermal disk will be 
essentially two-dimensional in the sense that vertical HE is maintained and the waves 
propagate as spiral density waves without any change in the disk scale height. It 
should therefore be noted that spiral waves in protoplanetary disks will only behave 
like pure density waves when the self-gravity of the disk is negligible and the EOS is 
nearly isothermal. 

To predict the height of a shock bore, consider a classical hydraulic jump. The 
height of a classical hydraulic jump can be predicted using the Froude number F of 
the pre-jump region, which measures whether the flow is rapid, F > 1, or tranquil, 
F < 1 (e.g., Massey 1970). Assuming that the jump is non-dissipative, the height of 
the post-jump flow /i2 can be found by 




(5.21) 



where hi is the height of the fluid in the pre-jump region. In the limit that F ^ 1, 
^2/^1 — V^F. This classical jump result can be used as a model for understanding 
the maximum height a shock bore reaches during the post-shock vertical expansion. 
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First, consider the definition of the Proude number, 
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F = 



(5.22) 



where g is the acceleration of gravity and Ui represents the pre-jump flow in the frame 
of the jump. Second, consider a non-self-gravitating disk. In this limit, one may write 
g —D,^z and hi ~ c^/Q, where Q is the circular speed of the gas, Cg is the midplane 
sound speed, and hi is the scale height of the disk in the pre-shock region. Using 
these relations in place of the corresponding terms in equation (5.22) reveals that for 
a non-self-gravitating disk, F — > Al. Any relation for the ratio of the scale heights 
before and after the expansion should have a behavior similar to equation (5.21), 
since F and M. are closely related. 

To relate the scale heights before and after a shock bore, consider equation (5.13) 
to be a measure of the overpressurization in the post-shock gas. Prom this perspec- 
tive, Qz in equation (5.16) represents a force field capable of doing work on the gas. 
Furthermore, assume that the difference in the gravitational potentials before the 
expansion and at the peak of the expansion measures work done in the expansion. 
Such an approximation allows us to write 




(5.23) 



/l2 

h~i 




(5.24) 
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where the RHS of equation (5.23) is the potential difference of the gas before and after 
expansion and the LHS of equation (5.23) is a measure of the work that can be done 
by the post-shock overpressure. Note that I neglect self-gravity in this approximation 
because for many astrophysical situations q > 100. This relation has a behavior 
similar to equation (5.21); when 7W = 1, h2/hi = 1, and, when A4 ^ 1, h2/hi ~ Ai. 
This does not mean that every fluid element is expected to jump to the new height 
given by equation (5.24). Instead, the scale height, as a characteristic height of the 
disk, will be changed. For example, in equation (5.18), material near the midplane, 
where 8^2/ dz — > 0, will hardly be affected by a single jump, while high altitude gas 
will have the strongest response. 

Equation (5.24) can also be derived using mass and momentum flux arguments 
(e.g., MC98). Mass conservation requires that Si^i = XI2M2, and momentum conser- 
vation requires that Pi — P2=T,iUi{u2 — ui), for Pi = J pidz. Again, assume that 
self-gravity is neghgible and that the structure before and after the shock is homol- 
ogous when in equilibrium; then p = Apoh"^, P = ABpoh^, and h = CS/po, where 
A, B, and C are shape factors that depend on the details of the model and po is the 
midplane density. By using mass and momentum conservation. 




(5.25) 
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Figure 5.2 A cartoon depicting the gas flow in a shock bore in the frame of the spiral 
shock inside corotation. The gas in the pre-shock region flows into the spiral shock 
(A). The shock (B) causes the material to be out of vertical force balance and a rapid 
expansion results (C). Due to spiral streaming and the loss of pressure confinement, 
some of the gas will flow back over the spiral wave and break onto the disk in the 
pre-shock region at a radius inward from where it originated (D). 
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Using equation (5.3) and assuming BC ^ 1, equation (5.25) becomes equation (5.24). 
Depending on the EOS, the assumption that BC ~ 1 may a bit inaccurate, but 
should introduce an error no larger than about 10% when using equation (5.24). For 
example, in an isentropic gas BC — r(l + m)r(l/2 + m)/[r(3/2 + m)r(m)], where 
m = — 1), which yields BC = 0.833 for 7 = 5/3, and BC — ^ 1, as 7 — 1. 

5.2.2 A Shock Bore in a Disk 

In a disk, shock bores turn into large waves that break onto the disk's surface, 
resulting in flow that is considerably more complex than described in §5.2.1. Shock 
bores will not only create vertical undulations in the disk, but will drive fluid elements 
to large radial excursions from their circular orbits and result in stirring the nebula, 
possibly mixing it. For reasons I explain below, inside the corotation radius of the 
spiral shock, these waves flow back over the spiral shock and break onto the pre-shock 
flow. This behavior is conflrmed in the numerical simulations in §5.3. I conjecture 
that shock bores will also result in breaking waves that crash onto the pre-shock flow 
outside the corotation radius, but I do not compute simulations outside corotation in 
this study. 

For the following discussion, consider the point of view from inside the corotation 
radius of the spiral shock. The development of breaking waves can be understood 
by recognizing two effects: When the gas crosses the shock front, the shock-normal 
component of a fluid element's velocity will be diminished by a factor given by the 
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inverse of equation (5.3), while the tangential component will be preserved, allowing 
the flow to be supersonic after the shock. This leads to streaming along the spiral 
arms, as demonstrated in streamline simulations of fluid elements in spiral galaxies 
(Roberts et al. 1979). In addition, when the gas expands upward, the pressure con- 
finement in the direction normal to the shock front is lost, and the material expands 
horizontally, causing some gas to flow back over the top of the shock. As the jumping 
gas moves inward and out over the pre-shock flow, it no longer has pressure support 
from underneath and breaks back onto the disk. The resulting morphology is a spiral 
pattern moving through the disk with breaking surface waves propagating along the 
disk's surface with the same pattern speed as the spiral wave (see Fig. 5.2). This 
morphology is only expected for a simple shock bore, i.e., there are no additional 
waves and shocks except for what are produced by cleanly defined spiral and break- 
ing waves. In a real disk with competing spiral waves, the behavior can be much 
more complex. 

In the simple case, there should also be a radius at which shock bores provide 
the strongest corrugation in the disk's surface. Consider the fluid elements in the 
disk to be essentially on Keplerian orbits. In a simple approximation, a shock bore 
produces a vertical perturbation to the fluid element's orbit, putting it on an inclined 
orbit leaving and then returning to the midplane in about half a revolution. In the 
inner disk, the orbit period of a fluid element is much shorter than the pattern period 
of a spiral wave. As a result, after a fluid element encounters the first spiral shock. 



112 CHAPTERS. SHOCK BORES 

the pattern will have moved only a small fraction of a pattern period by the time 
the fluid element encounters another arm; all fluid elements end up elevated between 
shock passages. As one moves radially outward in the disk, the advancement of 
the spiral shock becomes important and shock bores develop into breaking waves. 
However, as one moves toward the corotation radius, the shocks become weak, and 
shock bores are suppressed. 

5.2.3 Initial Models and Perturbations 

To study shock bores numerically, I generate spiral shocks by applying nonax- 
isymmetric perturbations to two equilibrium axisymmetric disks. I use the Mejia and 
cooled moderate disks as initial conditions (see Chapter 4) . The outer radius of each 
disk is scaled to 6 AU, and the central star is set to 1 Mq. The resulting inner holes 
for the two models are about 0.3 AU and 0.6 AU in radius. To relate these disks to 
Solar Nebula models, suppose that E(r) became r~^-^ in the unmodeled disk outside 
6 AU (see Lissauer 1987). Then the total disk mass would be 0.82 Mq out to 40 AU 
for the Mejia disk and 0.21 M© for the cooled moderate disk. The reader should keep 
in mind that these disks are only meant to demonstrate the dynamics of shock bores. 

Since shock bores can become very complex, especially when multiple jumps occur 
near each other, I stimulate a single, two- armed spiral wave. This is done using two 
methods. The first method forces spiral waves in the Mejia disk by adding a cos 2(f 
potential perturbation $p centered near 5 AU with a radial FWHM of about 0.5 AU. 
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This potential perturbation has the form 



$p (r, (fi;t) ^ A cos (2 {(p - ftpt)) cos^ (■ 



'Tr\r -rp\AR~^) , 



(5.26) 



where rp — 5 AU defines the pattern rotational speed Q,p, which is assumed to be 
Keplerian, A is some scaling parameter, and AR = 1 AU. This perturbation creates a 
well-defined, two-armed spiral that reaches all the way down to the central hole (see 
Fig. 5.3) and does not change the center of mass. Since the perturbation is localized 
to 5 AU, the change in the potential far from 5 AU, which stimulates the waves, can 
be thought of as due to concentrations of mass at the potential minima. In this way, 
the cos 293 potential perturbation behaves like two stubby spiral arms that could be 
produced by gravitational instabilities (see, e.g.. Boss & Durisen 2005a,b). 

For the cooled moderate disk, a difi^erent cos 20 perturbation is applied. Because 
the outer portion of the modeled disk has low mass, the stubby spiral arms produced 
by equation (5.26) do not generate strong shocks. Instead, I place two 2.5 Mj coro- 
tating point masses at r = 5.2 AU in the midplane of the disk, separated azimuthally 
by TT radians. Again, this method keeps the center of mass of the system at the cen- 
ter and creates two well-defined spiral arms, but it does allow for additional, small 
arms to form. Although I am treating the perturbers as point masses, one should 
not interpret them strictly as protoplanets, because any such perturber would open a 
gap within a few orbits (Bate et al. 2003). Instead, envision these point masses to be 
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transient clumps, formed by gravitational instabilities in a surface density enhanced 
ring, perhaps a dead zone (see Chapter 7). 



5.2.4 Thermal Physics 

The behavior of shock bores depends on the EOS, as demonstrated in §5.2.1. To 
test the predictions of the shock bore theory, the Mejia disk is evolved two ways. 

(1) An energy equation that includes some cooling plus heating by bulk viscosity is 
used with an EOS p = (7 — l)e, where e is the internal energy density and 7 = 5/3. 

(2) An isothermal EOS is used, where p(r, </?, z) — p{r, </?, z)TQ{r) and Tq represents 
the initial axisymmetric midplane temperature. The cooled moderate disk is evolved 
only with an energy equation. When the energy equation is used, a gentle volumetric 
cooling rate with a constant cooling time is applied to the disk to partially balance 
shock heating. Without any cooling, the disk heats up in less than a pattern period 
and surface waves are suppressed. The cooling time for the Mejia disk is set to four 
pattern periods at the perturbation radius. A constant cooling time is also used for 
the cooled moderate disk to balance shock heating, but with a cooling time of six 
pattern periods because the shocks are weaker. 
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Figure 5.3 Logarithmic surface density colorscale spanning 4.5 orders of magnitude 
for the Meji'a disk. Flow in this diagram is counterclockwise. 
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5.2.5 Fluid Element Tracer 

To investigate the shock parameters in both disks, the post-analysis fluid element 
tracer (§2.5) is employed, with data roughly every l/400th of an orp. Since the data 
storage demands are extremely high to achieve this type of spatial and temporal res- 
olution, the disks are first evolved at a low resolution, 128 aziumthal zones, until 
spiral shocks develop. Once the shocks form a quasi-steady morphology in the mid- 
plane, which takes about one or two pattern rotations, I switch to high resolution, 
512 azimuthal zones, and continue the simulations for a little more than half a pat- 
tern period. The isothermal calculation is only done at low resolution, because the 
low resolution is sufficient to demonstrate the general behavior of shock bores under 
isothermal conditions. 

5.3 Results 

In the simulations with heating and cooling (Energy Eq. in Table 5.1), the disk 
perturbations result in shock bores along the spiral arms. The jumping material 
creates very large waves that extend to nearly twice the height of the disk, and the 
only steady features in the simulations are the presence of the shock fronts and the 
shock bores. The shock bores vary in strength, and the nonlinearity of the waves 
creates transient features. 
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Figure 5.4 A radial cross-section portraying the shock at r = 2.5 AU for the Meji'a. 
The thick hnes are density contours corresponding to 3.5(-12), 3.5(-ll), 3.5(-10), 
7.1(-10), and 1.4(-9) g cc^^. The gray, shaded regions indicate shock heating corre- 
sponding to 5.6(-9), 2.2(-8 ), 9.0(-8), and 3.6(-7) erg cc~^ s~^. The arrows show the 
velocity of the gas, with each component scaled to its appropriate axis. 
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Figure 5.5 Cylindrical cross-section at r = 2.5 AU for the same time shown in Figures 
5.3 and 5.4. The density and heating contours are the same as in Figure 5.4. The 
arrows represent the flow of the gas in the frame of the potential perturbation and 
the color of the arrows represent the radial flow of material through this cross-section 
(blue toward the reader or central star and red away). This cross-section is shown from 
the perspective of an observer at the central star, and (p is measured counterclockwise 
from the 3 o'clock position in Figure 5.3. 
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These bores and breaking waves created by spiral shocks in disks are fully three- 
dimensional; looking at only a single slice of a wave's morphology and its correspond- 
ing gas flow can be both confusing and misleading. In the following discussion, I 
present two 2D-cuts, r-z and (f-z at different disk positions, as presented in Boley 
& Durisen (2006). When considered together, they form a comprehensible picture 
of waves in the disk. In the slices, contours delineate density and shock heating by 
artificial viscosity, and arrows represent gas velocity vectors. It is important to re- 
member that one cannot trace fluid element trajectories by connecting gas velocity 
vectors in any one plot. 

5.3.1 Meji'a Disk 

The cross-sections for the energy equation simulation show a very clean spiral 
shock at r = 2.5 AU, so I will describe it first. As a point of reference. Figure 5.3 
is the midplane density grayscale of the disk at the same time that the following 
cross-sections are shown. Figure 5.4 shows the velocity vectors of the gas in an 
r-z cut, which corresponds to 8 o'clock in Figure 5.3. The arrows represent the 
velocity components of the gas, with each component scaled to its axis, the thick 
contours represent the density structure, and the light contours with gray fill represent 
shock heating due to artificial viscosity. Near the midplane, where fluid elements are 
affected weakly by the shock bores, fluid element trajectories are roughly arranged in 
streamlines, creating the oval distortion seen in Figure 5.3, with the apastron near the 
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spiral shock. This is indicated by the gas velocity vectors near 2.4 AU in Figure 5.4. 
For r > 2.4 AU, the gas at all altitudes has already passed through a spiral shock at 
a larger radius and is flowing radially inward. The downward moving, high-altitude 
material that seems to appear from no origin is passing through this cross-section 
and is part of a coherent flow. For r < 2.4 AU, the low altitude material {z < 0.3 
AU) represents the outwardly moving pre-shock flow, while the mid- to high-altitude 
disk material has already passed through a spiral shock and has now developed into 
a large breaking wave. The breaking wave produces high-altitude shocks over the 
pre-shock flow and a large vortical flow, which appears to have some effect on the 
pre-shock flow even at low-disk altitudes. The disk height reaches a maximum near 
the spiral shock and is slightly less than twice its unperturbed height. The height of 
the shock bore will be discussed quantitatively in §5.4. 

Figure 5.5 is a ip-z cross-section at r = 2.5 AU, where (p is measured counterclock- 
wise from the 3 o'clock position in Figure 5.3, with velocities shown in the reference 
frame of the spiral wave, and complements Figure 5.4. In this flgure, the arrows are 
colored to represent the radial flow of the gas as seen by an observer positioned at the 
star. Before the shock, the gas is moving radially outward, and after the shock, it is 
moving radially inward. In addition, the mid- to high- altitude gas starts its vertical 
expansion at that interface. The sharp wall of gas at 2.5 AU X(p — 6 AU is formed by 
gas that jumped at a larger radius and is now falling inward through this cross-section. 
As this material crashes back onto the disk, it creates high-altitude shocks that result 
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in the expansion of the high-altiude gas as the post-shock region transitions into the 
pre-shock region of the next arm between 2.5 AUx0 = 6 to 8 AU. The gas begins to 
move radially outward again due to the roughly elliptical orbits of the gas. The wall 
of material evolves with time, and can develop a tube similar to a breaking surface 
wave. The structure of shock bores is time-dependent and transient^. 

Figure 5.6 shows the region 0.6 AU inward at about 6:30 o'clock in Figure 5.3. 
The breaking wave drastically affects the high density gas flow. A clear shock front 
cannot be defined on the coarse 3D grid, but strong shocks are everywhere in the 
wave. A large vortical flow reaches down to the midplane at r = 1.8 AU; another 
vortical flow is present in the high-altitude gas at r = 2.0 AU. These circulation cells 
are clearly related to the shock bores and waves. In the corresponding cylindrical 
cross-section. Figure 5.7, the complexity of the shock structure is again highlighted, 
and it is difficult to distinguish the shock bore, which is at about 1.9 X(p AU = 9 AU, 
due to the strong disturbance by gas that jumped at a larger radius and is plunging 
into the disk between 1.9 x ip AU = 7 and 8 AU. 

On the other hand, the shock bore shown in Figure 5.8 for r = 3.1 AU at about 
9:30 o'clock in Figure 5.3, which is near the Lindblad resonance of the imposed m = 2 
perturbation, shows the development of a very large wave that has yet to break onto 
the disk's surface. The cylindrical cross-section for this location. Figure 5.9, indicates 

that the wave is associated with a very sharp wall, similar to the 2.5 AU morphology. 

^An animation showing the time-dependence of shock bores is made available at 
http://hydro.astro.indiana.edu/westworld. Download the file Shock Front under the Movies link. 
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Figure 5.6 The same as Figure 5.4, except cutting across the shock at r = 1.9 AU. 
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Figure 5.7 Same as Figure 5.5, except for r = 1.9 AU. The shock morphology is much 
more complex and the shock bore cannot easily be defined by a single shock. Note 
the material, which jumped at a larger radius, plunging down into the disk between 
7 and 9 AU. 



5.3. RESULTS 125 

These figures show that at all disk radii portrayed, the disk material at high 
altitudes is moving differently from the material near the midplane, shocks are present 
at all disk altitudes, and the shock bores lead to large breaking waves, which create 
extensive vortical flows in the disk. The effect that these bores and waves have on the 
disk surface is emphasized by the isodensity surface contour shown in Figure 5.10. 

The high-mass disk was also evolved using an isothermal equation of state for 
the same potential perturbation. Figure 5.11 shows a (f-z cross-section for a shock 
near r = 2.5 AU. Instead of expanding at the shock, the gas compresses slightly and 
the height of the disk remains about the same. I speculate that the small peak of 
downward-moving, low density material just before the spiral shock is a remnant of 
transient features created by suddenly switching the EOS of the disk and by forcing a 
strong disk perturbation, and it is not part of the shock bore. The vertical compression 
of the gas at the shock front is in agreement with equation (5.18) because Jj <1 for 
an isothermal shock as given by equation (5.20). In an isothermal gas disk, a spiral 
shock will behave like a spiral density wave if self-gravity is unimportant, or it will 
compress vertically along the shock front if self-gravity is important. 
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Figure 5.8 The same as Figure 5.4, except cutting across the shock at r = 3.1 AU. 



5.3. RESULTS 127 




10 

3,1 AU M ^) 



Figure 5.9 Same as Figure 5.5, except for r = 3.1 AU. Even though this cross-section is 
the closest to the corotation radius, a strong shock and shock bores are still observed. 
This is probably due to the large radial motions of the fluid elements. 
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Figure 5.10 Isodcnsity surface contour for p = 3.54(— 10) g cc~^ shown at the same 
time as Figures 5.4-5.9. The view is from above the disk looking from about 5 o'clock 
to 11 o'clock in Fig. 5.3. The surface contour routine is unable to show breaking 
waves clearly. 
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Figure 5.11 A cylindrical cross-section at r = 2.5 AU for the isothermal Meji'a disk 
simulation. The density contours are the same as in Figure 5.5. The shock front does 
not cause rapid expansion in the post-shock region, but instead, causes a compression. 
I speculate that the small peak of downward-moving, low density material just before 
the spiral shock is a remnant of transient features created by suddenly switching the 
EOS of the disk and by forcing a strong disk perturbation. Although this wave is not 
associated with the compression caused by the shock bore, it is probably responsible 
for the slight undulatory morphology in the post-shock region. 
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Figure 5.12 The same as Figure 5.4, except for the moderate disk cutting across the 
shock at r = 2.7 AU. 
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Figure 5.13 The same as Figure 5.5, except for the moderate disk at r = 2.7 AU. 
The shock is weaker but still strong enough to induce an shock bore. In addition, the 
peak at 11 AU, which is moving radially outward, is not obviously associated with 
the shock bore. The nonlinear dynamics of waves in these disks is complex; the shock 
bores are not necessarily the whole story. 
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5.3.2 Cooled Moderate Disk 

The cooled moderate disk portrays a different shock morphology from the Mejia 
disk. The shocks are weak, and therefore, the shock bores are weaker and are localized. 
The wave structure in the disk may be more closely related to /-modes, as described 
by L098. Figure 5.12 shows the r-z cross section for the shock at r = 2.7 AU, which 
is at the same relative position as the 2.5 AU cross-section of the Mejia disk. There 
is surface corrugation and a weak shock bore at the spiral shock, but a breaking wave 
does not form. Figure 5.13, the cylindrical cross-section, shows that a strong vertical 
wall formed by inward falling material is, for the most part, absent. Nevertheless, the 
gas does expand vertically after entering the shock. 

5.4 Discussion 

5.4.1 Shock Bores 

Each of the three simulations shows different shock and wave morphologies. In 
the Mejia disk, strong shocks are produced, which lead to the formation of shock 
bores over most of the spiral wave. The spiral waves in the isothermal calculation 
propagate principally as spiral density waves. The cooled moderate disk only shows 
a bore-like morphology near r = 2.7 AU, and the disk has weak spiral shocks. It 
appears that shock bores are the extreme nonlinear outcome when /-modes develop 
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into strong spiral shocks. To understand how well the shock bore theory describes 
this nonlinear regime, I use the fluid element tracer to quantify shock parameters. 

Figure 5.14 shows a sample thermal history for a fluid element starting at r = 2.5 
AU, z — 0.176 AU, and — 0° in the Mejia disk. The fluid element is placed in the 
disk when the shocks first become strong, at the same time that the grid is switched 
to high resolution. The fluid element is followed for half a pattern period. The shock 
is seen clearly in the density and temperature/pressure plots. It is also noticeable 
in the u profile, where u is the shock normal speed in the frame of the pattern. To 
estimate u, it is assumed that the pitch angle i of the spiral arms is constant, nearly 
20° for these calculations. With this assumption 



Although this will not precisely measure u, and therefore the Mach number, it will 
give a reasonable value for most of the shocks in the disk. A modification to this 
method is made for studying shock strengths in disks for which the spiral waves are 
uncontrolled (Chapter 7). Figure 5.15 indicates the Mach number as well as the 
trajectory of the fluid element corresponding to Figure 5.14. Even though the fluid 
element starts at r = 2.5 AU, it is transported outward to 3.1 AU (Fig. 5.8) before 
it encounters the shock. The Mach number for the shock in Figure 5.8 is around 2.6, 
which roughly agrees with the change in density. The resulting change in the disk 



u = Vr COS i + {fl — Qp)r sin i. 



(5.27) 
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height should therefore be about 1.7, in the non-self-gravitating limit, or 1.6 assuming 
a g = 10; the actual vertical jump is about 1.5. Similar calculations indicate that 
the Mach number for the spiral shock shown in Figure 5.4, the radius where the fluid 
element is originally launched, is about 2.7. The jump in disk height again appears 
to be about 1.5-1.6. 

Even though the jump morphology in the r = 1.9 AU cylindrical cross-section of 
the Mejia disk is unclear, the shock bore description of spiral shocks predicts the new 
height for the inner disk fairly accurately. The Mach number for the shocks at r = 
1.9 AU is also about 2.7. According to equation (5.24), the new disk height should 
also be about 1.6-1.7 times higher for g — > cxo than the original height. The initial 
height of the disk is about 0.3 AU at r = 1.9 AU and the height of the disk is about 
0.5 AU in Figure 5.7, which matches the prediction well. 

The simple shock bore picture may be complicated by other wave dynamics near 
the Lindblad resonance at r = 3.1 AU (see Fig. 5.8 and 5.9), but shock bore theory still 
seems to predict the behavior of the shocks well in this region. Along the spiral shock, 
a bore is clearly seen. The fluid element tracer indicates that the Mach number for the 
shock at mid-disk altitudes is near 2.6, as stated above. The Mach numbers probably 
remain high at these radii because the fluid elements are on more highly elliptical 
orbits than in the inner radii, and the radial component to u becomes important, see 
equation (5.27), for these small pitch-angle spirals. 
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Figure 5.14 An example of a thermal history for a fluid element starting at r = 2.5 
AU, z = 0.176 AU, and (/? = 0°. Top left: Temperature (solid) and pressure (dotted) 
histories. Top right: Density history. Bottom left: Disk altitude trajectory (solid) and 
Vz (dashed). Although the vertical velocity does become negative and the vertical 
altitude reaches a local minimum while it is going through the shock, it quickly 
expands with a very sharp change in the vertical velocity as soon as it is in the post- 
shock region. Roughly, the shock bore is between the hash marks. Bottom right: The 
shock normal velocity, with respect to the global spiral shock, in the frame of the 
shock. 
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Figure 5.15 To complement Figure 5.14, the fluid element trajectories are shown 
above. Clockwise, starting in the upper left: r vs. z, r vs. time, the projection of the 
motion onto the midplane in Cartesian coordinates, and Mach number vs. time based 
on the P, p and u information presented in Figure (5.14). The open arrow heads 
show where the shock begins. Roughly, the shock bore is between the hash marks. 
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For the cooled moderate disk, the fluid element tracer indicates that the spiral 
wave forms a distinct shock near r = 2.7 AU. Since the shock is weak, u is unreliable, 
so I put the measured density change into equation (5.3) to flnd a Mach number of 
about 1.7. Assuming self-gravity is negligible, I expect that the maximum height of 
the jumping material to be about 1.3 times the initial scale height, which is fairly 
consistent with Figure 5.12. 

In discussing the changes of height in the disk, I have used the change in the 
density contour morphology to estimate differences. However, this can be misleading. 
To give a better quantitative description, deflne a local disk scale height to be 



Figure 5.16 plots this scale height for all three simulations at a similar radius. The 
scale heights in the pre-shock region for the energy and isothermal calculations are 
similar, but the isothermal disk shows a decrease in the scale height after the shock, 
as predicted, while the energy equation disk shows an increase. The change in disk 
scale height for the energy equation is about 1.5, compared to 1.6 or 1.7 predicted by 
equation (5.24) for the shock at this radius. The cooled moderate disk has a more 
gradual change in scale height than the other calculations, because of its weak shocks. 
The scale height changes by a factor of 1.3, close to what is predicted by equation 




(5.28) 



(5.24). 
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Figure 5.16 The scale height of each disk vs. azimuth. The Mejia disk, both energy 
and isothermal, are plotted for r = 2.5 AU and the cooled moderate disk is plotted 
for r = 2.7 AU. The shocks all occur between about 2 to 4 radians. 
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Figure 5.17 The results of tracing 1000 randomly selected fluid elements within an 
annulus for half a pattern period. Top left: The initial radial and azimuthal fluid 
element positions. Top right: The position of the same elements after half a pattern 
period. Bottom left: The initial vertical distribution. Bottom right: The final vertical 
distribution. Note that even material near the midplane is starting to show signs of 
stirring. 
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5.4.2 Processing and Mixing 

Large-scale vortical flows are present in both disks, which could result in the 
mixing of disk material radially and vertically. In addition, the high-altitude disk 
material can be moving at a much greater speed, and even in a different direction, than 
the mid- and low-altitude material. This makes it is possible to have a slip surface 
or vortex sheet at a Mach intersection, i.e., an intersection of two parallel flowing 
streams with different Mach numbers (e.g., Massey 1970). This will provide additional 
turbulence in the disk on a small scale and could influence the evolution of solids in 
the disk, e.g., size-sorting of chondrules (e.g., Cuzzi et al. 2001). Unfortunately, this 
effect cannot be modeled in these simulations since the typical cell size is ~ 10^^ cm, 
which is too large to study turbulence on the Kolmogorov scale where chondrules 
would be size-sorted (Cuzzi 2004). Nevertheless, the energy contained in any shock 
bore, or even a strong wave, is enough to provide the necessary chondrule size-sorting 
turbulence in a disk (see Boley et al. 2005). 

The large-scale vertical and radial mixing/stirring produced by the spiral shocks 
and shock bores can be investigated by using the fluid element tracer. Figure 5.17 
shows the results of tracking 1000 fluid elements for half a pattern period. The fluid 
elements are randomly distributed in disk volume between r = 1.44 and 3.36 AU and 
between z — and 0.24 AU. For this analysis I cannot distinguish between mixing 
and stirring of nebular material, mainly due to the short integration. Figure 5.17 
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demonstrates, however, that shock bores do at least produce significant stirring over 
tenths of an AU in only about 5.5 years. 

Episodic mixing and stirring of the nebula could be driven by spiral shocks that 
form from repetitious clump or arm formation over several million years. This would 
also provide short-lived shocks that could process solids throughout the nebula (Boss 
& Durisen 2005a,b), and may generate disk turbulence (Boss 2004b; Boley et al. 2005). 
Waves and shock bores could represent the major mechanism by which solids are 
processed and packaged into their parent bodies in protoplanetary disks. I return to 
this idea in Chapter 7. 

5.4.3 Shock Bores and Convection 

Convection can activate in a disk when (3 > (87 — 4)/(7 — 1) for an opacity law 
K = KqT'^, where 7 is the ratio of specific heats (Lin & Papaloizou 1980; Ruden & 
Pollack 1991). If this condition is satisfied, the vertical entropy gradient is driven neg- 
ative, and convective instabilities can grow and convective cells can be established for 
a range of disk radii (e.g., Lin & Papaloizou 1980; Rafikov 2007). However, these con- 
vective cells can also be easily disrupted by additional dynamics. Boley et al. (2006) 
found that for the M2004 simulation (see below), convective cells were dispersed by 
the onset of spiral wave activity. In contrast, other researchers claim that convection 
not only is active in very dynamic disks, but it has stark consequences for disk cool- 
ing. In particular. Boss (2004a) claims that convection leads to rapid coohng and 
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disk fragmentation. Moreover, Boss (2002) finds that disk cooling is insensitive to 
metaUicity due to convective fluxes. In this section, I explore the possibility of con- 
vection in unstable disks, and demonstrate that sudden vertical motions in dynamic 
disks are most likely shock bores and not convection cells. 

I search for convection in several unstable disk simulations by looking for regions 
where upward and downward motions are commensurate with regions of negative 
entropy gradients and are not associated with strong shocks, as strong shocks will 
result in a shock bore. Figure 5.18 shows the Mejia disk simulation discussed in 
this section. Within some regions with a negative vertical entropy gradient, vortical 
flows are seen. However, shocks are typically found in these regions as well. Such an 
alignment suggests that the regions are dynamic, i.e., they are not representative of 
regions where convection occurs in a stellar interior sense, but are regions that are 
associated with shock bores and waves. However, one does not expect for convection 
to be self-sustained in such a disk because the cooling time is ad hoc. To provide 
further comment, I investigate simulations with realistic cooling. 
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Figure 5.18 Similar to Figure 5.4, but at a different location in the Mejia disk with 
additional contours denoting superadiabatic regions. The thin contour lines indicate 
where the vertical entropy gradient is negative. 
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Figure 5.19 Convcction-likc motions in an optically thick protoplanctary disk model 
(see text in §5.6). The heavy curve roughly indicates the disk's photosphere, and the 
arrows, which are scaled to each axis and to the midplane density for each column, 
indicate the momentum density. Typical Mach numbers for the gas range between a 
few hundredths to a few tenths. Convection-like eddies are present throughout most 
of the disk. However, cooling times remain long in this disk, because ultimately, the 
energy must be radiated away. 



5.4. DISCUSSION 



145 



Vertical Momentum/Pgy Midplane p/p^^ 

] 

■ mm 

3456789 10 
r(AU) 

Figure 5.20 The left panel shows the vertical momentum density through the r — 
plane at 2; = 0.4 AU, scaled to the azimuthally averaged midplane density pav For 
comparison, the right panel shows the midplane density, scaled to Pav The vertical 
motions inside 6 AU are associated with the development of a spiral wave, and so 
it is difficult to determine whether such motions are related to a convective flow or 
to wave dynamics, in particular, weak shock bores. However, the strong vertical 
motions centered on 6.5 AU and again at about 8 AU have no corresponding density 
enhancement. These are also the regions that have large supcradiabatic regions, and 
so those vertical motions appear to be bona fide convective flows. 
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Consider the simulation presented by Boley et al. (2006; see Chapter 6). The 
initial disk is the Mejia disk scaled to 40 AU with a primary mass of 0.5 Mq, and 
the disk was evolved with the M2004 radiation scheme (see §2.2) under an ideal gas 
law with 7 = 5/3. Extended superadiabatic regions develop during the symmetric 
phase of the disk's evolution. In these areas, the average Mach speed in the convec- 
tive eddies {vz/cg) ~ 0.06, where {vz/cg) — J p \vz\ /cg dV/ J p dV and the integrals 
are evaluated over the volume spanning between 15 and 25 AU in r. The value of 
{vz/cs) fluctuates between 0.01 and 0.15 when each annulus of grid cells is evaluated 
separately (Ar = 1/6 AU), and some of the convective eddies result in nontrivial com- 
pressional heating through artificial viscosity. To evaluate energy transport by con- 
vection, I estimate the convective flux Ff. — —l/2cppTvzi (dlnT/dz — VaddlnP/d^;) 
(Lin & Papaloizou 1980) through each cell in the volume described above. Here, the 
mixing length i — min {z, P/fllzp) and Cp is the specific heat at constant pressure. 
By dividing the total internal energy within the half disk between 15 and 25 AU in r 
by the total convective energy loss rate for that same region, I find that the convective 
cooling time is about 1 ORP and that it is comparable to the radiative cooling time. 
However, this method measures the convective flux based on a superadiabaticity es- 
timate according to the mixing length formahsm of Lin Sz Papaloizou (1980). For a 
second measurement, I calculate the energy carried by convective motions through a 
plane that cuts through the volume at 2; ~ 1 AU. The convective flux through the ith 
cell is Fc\i — pVzCpAT\i, where AT is the difference between the actual temperature 
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at the cell center and the azimuthally averaged temperature for that r and z. I find 
a convective cooling time of about 2 ORPs with this method. According to either 
estimate of the convective flux, both of which are crude and uncertain, convection is 
efficient at redistributing energy in the disk's interior. However, these high convec- 
tive fluxes may be due to a combination of the random perturbation to the initial 
constant vertical entropy proflle, to the seeding of superadiabatic regions at the in- 
terior/atmosphere interface, which is typical of the M2004 and C2006 schemes (see 
Chapter 3), and to the exclusion of irradiation, which tends to produce stabilizing 
temperature stratiflcation. Moreover, as discussed below, the energy ultimately must 
be radiated away at the photosphere. 

The convection cells that are established in the axisymmetric phase are disrupted 
once nonaxisymmetry sets in and spiral waves dominate the dynamics. Vertical mo- 
tions in disks with strong spiral shocks are related to shock bores, which create high- 
altitude shocks and help stabilize the disk against convection. 

Because shock bores deposit energy into the upper layers of the disk, it is reason- 
able to ask whether shock bores themselves provide convection-hke cooling. This does 
not seem likely. Shock bores effectively reduce the temperature in the post-shock re- 
gion as the jumping gas expands vertically. This removes thermal energy in the 
post-shock region and deposits that energy in the upper layers via waves. Therefore, 
I do expect shock bores and possibly other wave effects (e.g., see Lubow & Pringle 
1993, Lubow & Ogilvie 1998, and Ogilvie 1998) to enhance disk cooling. Although 
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energy may be effectively removed from the post-shock region and allowed to radiate 
away quickly, the interior of the disk can still only cool by radiative diffusion because 
of the positive vertical entropy gradient established by the shock bores. Prom this 
point of view, shock bores and other wave effects limit the efficacy of spiral waves in 
heating the disk, but do not enhance midplane cooling. Distinguishing between shock 
bores and convection is not a point of semantics. Shock bores are born of large-scale 
shocks in a disk, while thermal convection and the criterion for convective instability 
are usually described in the context of a disk in vertical hydrostatic equilibrium. 

One objection to the above analysis of the Boley et al. (2006) disk, and also the 
similar simulation with the BDNL scheme (Boley et al. 2007c), is that the disk is 
not highly optically thick, with the midplane Rosseland vertically integrated optical 
depth Tm ^ 100. To address the role of convection in a highly optically thick disk, 
with Tm > 1000 for a large region of the disk, I use the fiat-Q model without the 
dense ring for a numerical experiment. 

As a reminder, the disk is 10 AU in radius, approximately 0.1 Mq disk mass, 
and surrounds 1 Mq star. I reset j — lA everywhere and force the opacity law 
K — (T/150 K)^ cm^ g~^. This opacity power law, along with the chosen 7, ensures 
that the disk should be convective as discussed above. The model is moderately 
stable to GIs, with Q f=:i 2 for most radii. Figure 5.19 shows that convection appears 
to be active in this model. There are large pockets of negative entropy gradients at 
mid-disk altitudes. However, are these motions truly convection, or are they related 
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to other wave phenomena? Figure 5.20 compares the vertical mass flux through 
a plane at ^ = 0.4 AU with the the midplane density structure, each normalized 
by the azimuthally averaged density for every radius. The midplane density plot 
indicates that some of the vertical motions might be explained, in part, by wave 
dynamics driven by the weak spiral structure in the disk. However, near 6 AU where 
the convection- like motions are the strongest, there is no indication that spiral wave 
activity is correlated with the vertical motions. 

Crude measurements of energy transport by vertical gas motions, with Fc — 
pCpVz^T, indicate that vertical motions are as important as radiation in transporting 
energy vertically. Even though a large fraction of the energy in the low to middle re- 
gions of the disk can be transported by vertical gas motion, I find icooi^ ~ 1000 near 5 
AU. For comparison, fragmentation occurs when icooi^ < 12 for a 7 = 7/5 gas (Rice 
et al. 2005). The energy must ultimately be radiated away near the photosphere, 
and so the cooling times remain long because the cooling time at the photosphere 
regulates convection. This result is consistent with Rafikov's (2007) analytic predic- 
tions, with my numerical tests (see Chapter 3), the numerical simulations of Boley et 
al. (2006, 2007c), and the numerical simulations of Nelson et al. (2000) and Nelson 
(2000). I note again that Nelson et al. assume a vertically isentropic density profile 
when calculating the cooling times for their 2D SPH calculations, which is similar to 
assuming efficient convection. 

In the simulations by authors who claim to see fast cooling due to convection (see 
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Boss 2004a; Mayer et al. 2007), the vertical motions are likely shock bores, and the fast 
cooling due to "convection," I argue, is likely a result of poor boundary conditions 
for radiation physics or sudden changes in numerical parameters (see Chapter 7). 
As discussed by Nelson (2006), Boley et al. (2007c), and here, proper treatment 
of radiation physics, especially near the photosphere, is crucial for estimating proper 
cooling times. In fact, in order to lower the coohng times for the flat-Q model decribed 
above to those expected to lead to fragmentation, the effective temperature would 
need to be approximately equal to the midplane temperature, which is about three 
times the actual effective temperature at r = 5 AU. 

5.5 Conclusions 

5.5.1 Shock Bores And Waves 

When a strong spiral shock develops in a disk, the shock's highly nonhnear be- 
havior creates a shock bore, where the loss of vertical force balance in the post-shock 
region results in a rapid vertical expansion of the gas. The resulting shock structure 
covers a large range of disk altitudes. The shocks are not limited to the main spi- 
ral wave itself. Jumping gas can fall back onto the pre-spiral shock gas producing 
breaking waves, vortical flows, and mid- and high-altitude shocks. The analytic ap- 
proximations in §5.3 do an accurate job of predicting the height of the disk in the 
post-shock region as measured in 3D hydrodynamics simulations. 
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Protoplanetary disks with clumps and spiral arms can have very complex wave 
dynamics. As described by several authors (e.g., Bate et al. 2003), waves can transport 
angular momentum and lead to the formation of gaps. In addition, as suggested by 
these simulations, waves may provide turbulence in the disk as well as large vortical 
flows, some of which may extend down to the midplane. The combination of radial 
transport and vortical flows should work to stir, if not mix, the gas and the lighter 
solids over tenths of an AU in only an orbit period. Moreover, the wave dynamics 
should lead to copious shocks in the disk. 

I find that convection in disks is not a viable mechanism for significantly reducing 
the cooling times to those required for fragmentation. Moreover, vertical motions 
associated with shocks are hkely shock bores and not convective flows. 

5.5.2 Bearing On Reality 

As demonstrated here, in Boley et al. (2005), and in Boley & Durisen (2006), a 
perturbation must be strong to produce dynamic waves and shock bores in proto- 
planetary disks. However, such strong perturbations may be plausible. The onset of 
gravitational instabilities could create massive spiral arms and/or clumps that could 
drive strong spiral shocks throughout the disk, including the inner regions (Boss & 
Durisen 2005a,b). Even if the disk is gravitationally stable as a whole, there still may 
be regions in the disk, e.g., a dense ring or annulus formed by a dead zone (Gam- 
mie 1996), where gravitational instabilities could set in, produce a burst of spiral 
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wave activity, heat the region to stabihty, and possibly later repeat the process (e.g., 
Armitage et al. 2001). Such a transient mechanism could provide strong enough 
perturbations to drive strong waves and shock bores without opening a gap. The 
simulations presented here are not self-consistent inasmuch as a sudden perturbation 
is introduced and the simulations are not evolved for a long stretch of time. These 
simulations do, however, demonstrate that shock bores occur in strong spiral waves, 
and they provide a basis for comparison with much more complex simulations of GI 
active disks. 



Chapter 6 



GRAVITATIONAL DISK 
INSTABILITY 

In this section, I describe general results of two disk instability studies, namely the 
simulations presented in Boley et al. (2006) and Boley et al. (2007c). I find that the 
GI activity in these disks relaxes to a thermally self-regulating state, where heating 
balances radiative cooling. In this state, the effective Shakura & Sunyaev (1973) 
a ~ 10~^ due to Gl-driven torques, the transport is dominated by global modes, and 
the details of the mass transport and GI strength-mode spectrum is sensitive to the 
treatment of radiation physics. The cooling rates are too long for disk fragmentation 
to occur, and none is observed. 

The Boley et al. (2006; hereafter B2006) calculation is the same as the Mejfa 
(2004; hereafter M2004) radiative transfer simulation with no irradiation, which was 
evolved with the M2004 radiation scheme (§2.4.1; see also Mejfa 2004, §4.1.2). For 
the Boley et al. (2007c; hereafter BDNL) simulation, the M2004 disk was restarted 
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just after the burst, but evolved with the BDNL radiation transfer scheme (§2.4.2) in 
the standard version (SV) of the code (§2.1) to study how sensitive the evolution is 
to the treatment of radiation physics. 

6.1 Initial Conditions and Method 

The initial model is the Mejia disk (see Chapter 4), with the disk radius scaled 
so that it extends from 2.3 to 40.0 AU. The star is scaled to 0.5 Mq, and so the 
disk mass is 0.07 Mq. For these simulations, 1 orp (outer rotation period) 253 yr, 
and represents the initial gas orbital period at r = 33 AU. The M2004 simulation 
was run from the initial axisymmetric Meji'a disk, while the BDNL simulation was 
restarted from the M2004 simulation at about 6.5 orp. Restarting just after the 
burst was a compromise between lowering the numerical cost of the simulation and 
avoiding transient features introduced by the switch in radiation physics routines as 
the disk approaches an asymptotic state. For both simulations, due to an error with 
the inclusion of helium in the solar mix, the mean molecular weight used for the 
temperature and pressure ranges in the simulation is around 2.7 when it should be 
close to 2.3. This introduces a systematic offset no larger than about 16% into the 
temperature, which directly affects the cooling rates. However, the opacity law is 
roughly quadratic in temperature, which means the flux calculated from the diffusion 
approximation [F ~ T^^/i^) is roughly quadratic in temperature, too. The error in 
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the mean molecular weight should be an error that typically enhances the cooling 
by 1.16^ or about 4/3. Because I expect the cooling to be enhanced, fragmentation 
should be more likely in the M2004 and BDNL simulations than it would be with 
the correct mean molecular weight. I also note that for the BDNL simulation, only 
Rosseland opacities were used because the interpolation described in §2.4.4 was still 
being developed at the start of the BDNL simulation. 

The temperatures in these simulations never reach the dust sublimation temper- 
ature, near T ~ 1400 K (MuzeroUe et al. 2003), so the opacity is mostly due to 
dust. D'Alessio (2001) opacities are used, with minimum and maximum grain sizes 
Omin = 0.005/xm and Omax = lAini, respectively (see Mejia 2004 or Appendix A of 
B2006) and an interstellar grain size power law dn — a~^^^da. The minimum allowed 
temperature at all times is 3 K to simulate radiating into empty space. However, 
it should be noted that there is an inconsistency regarding the treatment of the hy- 
drodynamics and the radiation physics. In the SV, wherever the temperature drops 
below 3 K, the temperature is reset without correcting the specific internal energy, 
leading to an inconsistent pressure. This inconsistency is corrected in CHYMERA by 
resetting the specific internal energy of the cell as well as the temperature. 

Both simulations are evolved approximately up to 4000 yr, i.e., about 16 orp. 
The resolution used is the same as the constant cooling time simulation presented 
in Mejia et al. (2005). For the M2004 simulation, the z direction is fixed at 32 
zones. As indicated by the radiation transfer tests (Chapter 3), the BDNL scheme 
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leads to a more extended atmosphere than the M2004 scheme, and so the vertical 
direction was extended to 64 zones. The starting model fills the initial grid, namely 
(r, (f), z) = (256, 128, 32), and so when the disk rapidly expands during the burst, the 
grid is extended to 512 in r, while keeping the same cell resolution of Ar = Az — 1/6 
AU. M2004 and BDNL are also run on a high azimuthal resolution grid, (r, (p, z) = 
(512, 512, 32 or 64) to test for fragmentation (see §5.2). A random cell-to-cell density 
perturbation of amplitude \Ap/p\ = 10~'^ is applied at the very first step of the M2004 
run, which allows spiral modes to grow from the background noise as the disk cools. 

6.2 Disk Structures 

Figure 6.1 shows surface density images of the M2004 simulation. The disk under- 
goes the same four evolutionary phases described for the disks in Pickett et al. (2003), 
Mejia (2004), Mejia et al. (2005), and Cai (2006), namely the axisymmetric, burst, 
adjustment, and asymptotic phases. The axisymmetric phase is the initial cooling 
phase, before the instabilties activate. The burst phase describes the often violent 
onset of GIs. The adjustment or transition phase is the relaxation to the thermally 
self-regulating Gl-active phase during which heating is balanced by radiative cooling, 
i.e., the asymptotic phase. 

During the axisymmetric phase, the disk begins to pulsate radially. The pulsation 
is probably due to a combination of the loss of radial pressure support due to cooling 
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and of the absence of a mechanism for redistributing angular momentum. A dense, 
thin ring, in which Q < 1, forms in the outer disk just before the burst. Once 
nonaxisymmetry develops, the ring becomes part of the spiral structure. Between 2 
and 3 orp, GIs fully activate, and the disk rapidly expands as open spiral structure 
efficiently transports angular momentum outward. During the adjustment phase, the 
disk oscillates in size to a much greater extent than noted in the axisymmetric phase. 
In addition, the power in the now highly non-linear GI activity begins to dampen in 
response to shock heating and work done by gravity, and the disk transitions to a 
thermally self-regulating state. 

One-armed structure does develop during the adjustment phase, but because the 
star is held fixed at the origin, some of the dynamics are not accurately treated. 
Methods like moving the star to a location that brings the center of mass back to the 
center of the grid (Boss 1998) are not employed because this might incorrectly treat 
the dynamics as well. A preliminary leap-frog method has recently been developed 
to release the star for the Wengen test 4 simulations (Mayer et al. 2007, in prepa- 
ration), and a more accurate method is being developed for standard incorporation 
into CHYMERA, but the simulations presented in M2004, B2006, and BDNL were 
completed before such routines were developed. 
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Figure 6.1 Evolution of the M2004 simulation. All images show surface densities 
in logarithmic colorscale spanning 4.5 orders of magnitude, with the abscissa and 
ordinate in AU. The computational grid is larger than what is shown in the image 
(80 AU), and significant amounts of mass are not leaving the grid. A movie of the 
simulation is available at http://hydro.astro.indiana.edu/westworld under the Movie 
s link. 
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Figure 6.2 Comparison between the evolution of the M2004 (middle row) and BDNL 
(bottom row) simulations. The scaling is the same as in Figure 6.1. Although fairly 
modest, the differences between the simulations are discernible; the spiral structure 
in the BDNL disk is more washed out than in the M2004 disk. 
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;ure 6.3 Comparison between the evolution of the M2004 (middle row) and BDNL 



(bottom row) simulations: merdional volume density slices at 3 o'clock in Figure 6.2. 
The images show logarithmic density on a colorscale spanning 4.5 orders of magnitude. 
The major differences between the images are that the BDNL disk is more collapsed 
in the inner disk than the M2004 disk and is more flared than the M2004 everywhere 
else. 
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Figure 6.4 Fourier amplitude spectrum for the M2004 disk (black) and the BDNL 
disk (red), time averaged over the last two orp of each simulation. The bars represent 
typical fluctuations over the two-orp period. The M2004 profile has larger amplitudes 
everywhere except for m = 2, which is consistent with the M2004 disk being more non- 
axisymmetric. Both spectra can be fit with a functional form of ~ {m'^ + mQ)"", 
where n ^ 1.5. This behavior at large m may be indicative of gravitoturbulence, 
nonlinear mode coupling, or both. 
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Although not extreme, quahtative differences between the M2004 and BDNL sim- 
ulations are discernible in the surface density plots (Fig. 6.2) and meridional slices 
(f) Ki (Fig. 6.3). The M2004 simulation has more pronounced spiral structure than 
the BDNL simulation throughout most of the disk. Also, the BDNL simulation is 
more extended in the vertical direction because the disk's atmosphere is hotter, as 
expected from the tests. The ring that forms at the innermost edge of the BDNL 
simulation is due to poor numerical resolution; the mass at the edge cooled until 
contained within one cell. The M2004 simulation does not exhibit similar behavior 
because the scheme is slighly more dependent on resolution (Chapter 3). 

In order to quantify the structural differences, I compute the global Fourier ampli- 
tude spectrum for each simulation, where the sum over the amplitudes is a measure 
of the nonaxisymmetric structure in the disk and the spectrum is indicative of the 
dominate modes in the disk. I compute the time-averaged Fourier component (Am) 
for m-arm structure by 

^ J pordrdz ' 

where po is the axisymmetic density component and pm is the total Fourier amplitude 
of the cos(m0) and sin(m0) density component. The time-average is calculated by 
finding A^ for a large number of snapshots over the last two orps. The summed global 
Fourier amphtude = X)m=2 i^m) = 1-4 for the M2004 disk, while {A+) = 1.1 
for the BDNL disk, which indicates that the M2004 disk is more nonaxisymmetric. 
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I exclude m — 1 from the summation because the star is kept fixed. The difference 
between the sums is also depicted by the Fourier spectrum (Fig. 6.4). The M2004 
disk has larger amplitudes everywhere except for m = 2, which is consistent with 
the qualitative differences portrayed in Figures 6.2 and 6.3. Including long-range 
transport of radiation in the BDNL scheme has resulted in weaker GIs. 

As described in B2006, a curve with the functional form A^^ ~ (m^ + ml)""' can 
be fit to the data, where n ~ 1.6 and mo ~ 7.5. The BDNL disk is also consistent 
with this functional form, and both disks roughly follow n 1.5. The similar slopes 
at large m may be indicative of gravitoturbulence (Gammie 2001) or nonlinear mode 
coupling (e.g., Laughlin et al. 1998). The meaning of the power spectrum remains 
unsatisfactorily understood. 

6.3 Disk Energetics 

Figure 6.5 shows the evolution of the internal energy for each disk. There is a 
precipitous drop in energy over about an orp (253 yr) after switching the radiation 
transport schemes. However, the schemes roughly follow each other after the drop, 
with the BDNL profile having a more shallow slope than the M2004 profile. The 
effective temperature profiles (Fig. 6.6) are also similar; each profile can be fit by 
an exponential. As discussed in B2006, this profile is very steep when compared 
with observationally derived profiles that are based on spectral energy distribution 
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measurements. Observational studies typically find that for Tgff ~ r~*, q 0.4-0.8 
(e.g., Beckwith et al. 1990, Kitamura et al. 2002, and DuUemond et al. 2007; see also 
theoretical models by Miyake & Nakagawa 1995). The exponential profile is hkely a 
result of the exclusion of stellar irradiation in these simulations. 

For each simulation, I calculate the time-averaged cooling time 
tcooi — J ^dV/ / I V ■ F I dV for each annulus on the grid, where e is the inter- 
nal energy density of the gas and V • F is the radiative cooling. The temporal average 
is taken to be about the last six orp of evolution for the M2004 simulation and about 
the last 5 orp for the BDNL simulation. In Figure 6.7, 1 compare ^cooi^ curves for each 
disk, where fl is the angular speed of the gas. The cooling time is much longer for 
r < 35 AU in the BDNL disk than it is for the M2004 disk. This is hkely due to a com- 
bination of the different opacities used by the two routines and of the free-streaming 
approximation, employed by the M2004 algorithm, in regions where long-range ra- 
diation couphng matters (see Chapter 3). Outside r ~ 35 AU, the curves converge. 
The longer cooling times are consistent with the washed out structure in the BDNL 
simulation. For both disks, the cooling times are well above the fragmentation crite- 
rion icooi^ < 3 to 6 for a 7 = 5/3 gas (Gammie 2001; Rice et al. 2003), so I expect 
neither disk to fragment. Regardless, I ran the BDNL simulation at 512 azimuthal 
divisions between 10 and 11 orp (2530 and 2783 yr) to test for fragmentation, when 
the total internal energy reaches its minimum, and found no signs of fragmentation, 
as one expects from the long cooling times. 
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Figure 6.5 Internal energy for BDNL (heavy curve) and M2004 (light curve), normal- 
ized by the M2004 initial value. The precipitous drop at 7 orp is a result of suddenly 
switching radiation schemes. Between about 7 and 10.5 orp (1671 and 2657 yr) the 
curves approximately track each other. 
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Figure 6.6 Effective temperature profiles for the BDNL (heavy curve, time-averaged 
over about the last 5 orp or 1165 yr and M2004 simulations (light curve, time-averaged 
over the last 6 orp or 1418 yr). Both follow an exponential profile, and are reasonably 
consistent. Their departure from a rough r~^/^ effective temperature profile is likely 
due to exclusion of stellar irradiation. 
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Figure 6.7 Cooling time curves scaled by the local angular speed for the BDNL (heavy 
curve) and M2004 (light curve) simulations for the time-averaged periods over about 
the last 5 and 6 orp, respectively. Both curves are relatively consistent for r > 35 AU, 
but they depart for inner radii. This is likely due to a combination of the different 
opacities used by the two routines and of the free-streaming approximation, employed 
by the M2004 algorithm, in regions where long-range radiation coupling matters. 
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6.4 Angular Momentum Transport 

In this section, I compare the angular momentum transport in each disk by analyz- 
ing the gravitational torque on the inner disk due to the outer disk and by measuring 
the effective Shakura & Sunyaev (1973) a. As discussed in §2.5.2, the gravitational 
torque is calculated by 



where $ is the gravitational potential, x is the position vector, and the integral is 
over the volume contained inside r = R. For this analysis, I am concerned with the 
vertical component C^. The time-averaged torque, averaged over the last six orp for 
the M2004 disk and about the last 5 orp for BDNL disk, is shown for each simulation 
in Figure 6.8. The solid curves represent the torque profiles, with the heavy curve 
indicating the BDNL disk. The dashed curves show the mass flux for each disk with 
arbitrary but consistent scaling; the peak mass flux M — fewxlO"^ Mq yr~^ for each 
disk. The torques are of the same magnitude, but the torque profiles are noticeably 
different. Based on the {Am) plots and the visual differences in disk structures, 
the M2004 disk has a more complex morphology and stronger modes. The multiple 
peaks in the M2004 torque profile are another indication of this complex morphology 
and competing global, dominant modes. The BDNL torque profile also has multiple 
extrema, but the variations are not as extreme. The hydrodynamical stress (Chapter 
2) may play an important role in the BDNL disk. 




V{R) 
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Figure 6.8 The negative of the gravitational torque (sohd curves) and mass flux 
(dashed curves) profiles for BDNL (heavy curves), time-averaged over the last 5 orp, 
and M2004 (light curves), time-averaged over the last 6 orp. BDNL's torque profile 
shows one strong peak and several minor peaks, while the M2004 torque profile has 
two very strong peaks. The mass fluxes for each disk (arbitrary but consistent scaling) 
are consistent in magnitude. 
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Figure 6.9 Effective a profiles for BDNL (fieavy curves) and M2004 (ligfit curves). 
Tfie solid curves indicate the effective a derived from the torque profiles, and the 
dashed curves indicate the predicted a based on an a disk prescription, derived from 
the tcooi^ profiles (Gammie 2001) in Figure 6.7 with the assumption of negligible 
self-gravity (see text). Both disks roughly follow the predicted a over a large range 
of radii. 
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Figure 6.10 Same as Figure 6.9, but the BDNL simulation has only vertical radiation 
transport. 
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Both disks have comphcated mass flux proflles. The principal inflow/outflow 
boundary in the M2004 simulation is at r ~ 26 AU. The BDNL disk, by contrast, 
has two main inflow/outflow boundaries. The r ~ 15 AU boundary corresponds to 
the peak torque in the BDNL disk, and I refer to this as the principal inflow/outflow 
boundary because the mass fluxes are the highest near it. Roughly, the peak in 
each torque proflle aligns with the principal inflow/outflow boundary. The agreement 
is imprecise, because the mass flux average is based on differencing mass cylinders 
at different times, which yields a time average based on the second-order mass flux 
integration. The torques are derived in post-analysis calculations, and so the temporal 
sampling is much sparser. Moreover, the mass fluxes are highly variable with time, 
and averages over slightly different time periods can result in different mass flux 
proflles. However, major inflow/outflow transitions are usually near torque proflle 
extrema. The mass fluxes for r > 40 AU are complicated by pulsations that begin 
just before the disk bursts and continue throughout the evolution. 

I calculate an effective a by relating the vertically integrated gravitational stress 
% to the gravitational torque through % — — C^/27rr^ (equation [2.42]). An effective 
a is then given by (Gammie 2001) 



where the brackets indicate an azimuthally averaged quantity, c is the adiabatic mid- 



dlnQ -1 % 



(6.3) 



a = 



dlnr (c^5]) ' 
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plane sound speed, and E is the surface density. The comparison between the effective 
a for the BDNL and M2004 disks is shown in Figure 6.9. The profiles are of similar 
magnitude everywhere, but with BDNL being significantly lower between about 20 
and 36 AU. 

I also show in Figure 6.9 the a one would expect for an a disk (see Gammie 2001; 
equation [2.41]) 



where 72D is the two-dimensional adiabatic index. For 7 = 5/3, 720 ~ 1.8 in a 
strongly self-gravitating disk and 1.5 in a non-self-gravitating disk (Gammie 2001). 
I use the tcooi^ profiles (Fig. 6.7) in equation (6.4) to plot the anticipated a for a 
local model in the non-self-gravitating limit. For most radii, both disks are roughly 
consistent with this a prescription, and a is roughly constant between 20 and 35 
AU. The main difference between the two simulations is the lower a in the BDNL 
simulation, which is consistent with the longer cooling times. Recall that the free- 
streaming approximation is used in the M2004 scheme wherever r < 2/3. This can 
lead to additional coohng if r > a few hundredths (see Fig. 3.4); the use of only the 
Rosseland mean opacity in the BDNL simulation likely plays a role as well. 

Although the overall evolutions are in rough agreement, they demonstrate sen- 
sitivity to the details of radiation transport. By including the long-range effects of 
radiative transfer and by using different opacities for different optical depth regimes. 




(6.4) 
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the BDNL disk shows less structure, is more flared, and has a lower effective a that 
deviates slightly more from equation (6.4). These differences demonstrate the need 
for a radiation algorithm that includes the long-range effects of radiative transfer in 
all three-dimensions, which will be missed by diffusion approximations and is missed 
in the BDNL scheme in the r and directions. To illustrate the importance of ra- 
diative transport in all three directions, I show in Figure 6.10 the effective a proflle 
for a disk that was evolved with the BDNL radiative routine, but with the radial and 
azimuthal diffusive radiation transport turned off. According to this plot, without 
any r and transport, I would surmise that the effective a deviates strongly from 
the predicted a. 

6.5 Mass Redistribution 

Once GIs active, mass in the M2004 disk is quickly reordered. During the burst, 
roughly between 2 and 4 orp (500-2000 yr), accretion rates are as high as a 
fewxlO"^ Mq yr~^. During the adjustment phase (6-10 orp or 1500-2500 yr), the 
accretion rates decrease by an order of magnitude, and for the asymptotic phase 
(10-16 orp or 2500-4000 yr), the typical accretion rate decreases by another order 
of magnitude. The mass inflow/outflow boundary varies with time, and multiple 
inflow/outflow boundaries develop. 

Figure 6.11 (bottom) also illustrates how the r^^-^ surface density profile of the ini- 
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tial disk is lost, and approximately relaxes to a Gaussian of the form E = EqIO"^^/''''^^ 
for r > 20 AU (dashed curve). The least-squares fit in log-linear space with the 
dummy variable x — r'^ yields Tg = 46.7 AU when fit between r — [20, 60] AU. How- 
ever, on the intervals r = [20, 43] AU and r = [43, 60] AU the surface density profile 
seems to follow two different power laws E ~ r^*^. For the inner interval, = —1.93, 
and for the outer interval, u — —5.97. The Spearman correlation coefficients are 
R = —0.992, —0.986, —0.985 for the exponential, the inner power law, and the outer 
power law, respectively. 

The disk inside 20 AU is not well-described by any simple function. There is a 
ring at 7.5 AU, and another seems to be forming at 10.5 AU, both resembling those 
that appeared in the constant cooling simulation with the Mejia disk (Mejia 2004; 
Mejfa et al. 2005). Although the exact mechanism is unclear, these rings appear to be 
a numerical artifact that results from a combination between poor vertical resolution 
and keeping the star fixed. However, there is also a broad, radial mass concentration 
at about 15 AU due to tightly wrapped spiral arms, which is likely a real feature. I 
therefore speculate that if there are regions in disks that are kept hot for physical 
reasons, ring-like mass concentrations could grow. 
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Figure 6.11 Global mass redistribution for the M2004 disk. The top panel shows the 
average mass transport rate calculated by differences in the total mass fraction as a 
function of radius for three different times: between 2 and 4 orp, 4 and 10 orp, and 
10 and 16 orp. The red curve is scaled to the left ordinate while the blue and green 
curves are scaled to the right ordinate. The bottom panel shows the surface density 
as a function of radius for the initial disk (dotted) and the final state (solid), which 
reflects the density profile during the asymptotic state. The dashed curve shows that 
the density profile for r > 20 AU follows a Gaussian distribution as described in the 
text. However, the surface density profile can also be broken down into two power 
laws (blue curves). (Same as in B2006.) 
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Figure 6.12 Comparison between the surface density profiles of the M2004 and BDNL 
disks. Each profile is consistent with a Gaussian profile for r > 20 AU. 
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Figure 6.13 Toomre Q{r) for the M2004 and BDNL simulations. Both profiles are for 
snapshots at the end of the simulation. However, each profile evolves little between 
about 10 orp and 16 orp. The fine indicates Q — 1.5. 
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The BDNL simulation's mass distribution is very similar to that of the M2004 
simulation, despite the difference in the strength of GIs in each disk. Because the 
BDNL disk was restarted after the burst and the burst is responsible for reordering 
mass quickly, it is expected for the mass profiles to be fairly similar. As displayed in 
Figure 6.12, both surface density profiles follow Gaussians over a large radial extent. 

6.6 Thermal Balance 

The Toomre Q evolves httle between about 10 orp (2530 yr) and the end of the 
simulations (Fig. 6.13). As noted by M2004 and B2006, the average Q in the M2004 
disk between 14 and 43 AU is 1.5, while the average Q between the same region is 

1.7 for the BDNL disk. The BDNL disk having a higher Q than the M2004 disk 
is consistent with the previously discussed assessments. In addition, the steady Q 
indicates that both disks are in an asymptotic-like state, as was noted for the Mejia 
et al. (2005) constant cooling time simulation, and so the GIs in these disks have led 
to a thermally self-regulating state for which radiative cooling is balanced by work 
and shock dissipation. The distinction asymptotic- ZzA;e is made because the cooling 
times in both disks are continually adjusting, and so the behavior of GIs will evolve 
with the disk structure. To estimate the lifetime of the asymptotic phases for each 
disk, I note that the energy loss for the last few orp of evolution in both disks is 
roughly constant. Based on a linear fit to the energy curves for t > 14 orp (3540 yr; 
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Fig. 6.5), the M2004 disk can sustain the asymptotic phase for about 18000 yr, while 
the BDNL disk can sustain the same phase for about 25300 yr. For both disks, this 
asymptotic phase is relatively short-lived compared with typical total disk lifetimes 
(e.g., Hartmann 2005). However, effects hke grain growth and irradiation would likely 
alter these calculations significantly. The simulations presented by Cai et al. (2006, 
2007) indicate that grain growth and irradiation weaken GIs, which could result in 
slower evolution. Moreover, a relatively small disk (radially) is simulated in these 
calcuations, and much more extended disks may have different evolution timescales. 

6.7 Coherent Structure 

A periodogram analysis (Horne & Baliunas 1986) can be used to look for period- 
icities in data sets. When applied to the cosine of the phase angles of the radially 
dependent Fourier components of the azimuthal density distribution, the periodogram 
is a powerful tool for finding persistent, coherent structure in disks. In this case, the 
structure corresponds to m-arm spirals with steady pattern speeds. Structure that 
maintains the same phase for an extended period of time will have more power than 
noise or short-lived phase coherent structures. I should note that power in a mode 
on the periodogram does not necessarily mean that the mode is particularly strong, 
in the (^m) sense, although the two are often correlated. It means rather that the 
pattern is highly coherent. 
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Figure 6.14 shows periodograms for Fourier components m — 1, 2, 3, and 4 in the 
M2004 disk. Overplotted are the corotation speed, derived from the mass-weighted 
orbital motion of the gas, flw{r), and the inner and outer Lindblad resonances (ILR 
and OLR), which are calculated from the epicyclic frequency k — (2Q^dQ^/dr + 
40^)^/^. There is significant power in mode m = 1 at various points at corotation, and 
at the OLR at about 2 orp~^. Because the power is fairly locahzed at the resonances, 
the m — 1 power is likely a result of a combination between a superposition of modes 
(see below) and structure that arises from keeping the primary fixed. For 171 — 2, 
there are three strong swathes of power, with noise, corresponding to 0.8, 1.5, and 
2 orp~^. The corotation radius for the 1.5-orp^^ swath is commensurate with the 
ILR of the 0.8-orp~^ swath and with the OLR of the 2-orp~^ swath. In addition, 
these swathes of power are coincident with mass flux and torque profile features in 
Figure 6.8. The inflow/outflow boundary at r ~ 26 AU ahgns with corotation of the 
1.5-orp~^ swath, which also aligns with the outer peak in the torque profile. The 
2-orp~^ swath aligns with the dip in the mass flux proflle and with the inner peak in 
the torque proflle. 

In the 171 — 3 periodogram, the coherent structure is distinct. There are three 
strong signals at 0.8, 1.5, and 2 orp^^, as in the m — 2 periodogram, and there is 
additional structure at about 4.4 orp~^. For m — A, almost all coherent structure is 
washed out. Finally, I note that the m = 1 power along corotation is commensurate 
with the major m — 2 and 3 swathes of power. This is likely an indication of mode 
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coupling between the low-order modes in the disk (e.g., Laughhn & Korchigan 1996; 
Laughhn et al. 1998). 

Figure 6.15 is similar to Figure 6.14, but for the BDNL disk. Again, there is power 
in the m — 1 Fourier component along corotation. However, there is fairly coherent 
structure at 2 orp^^ between 15 and 40 AU, with a break of about 5 AU. This will be 
discussed more in a moment. For m = 2, there is very persistent structure between 
10 and 30 AU, with a well-defined pattern at about 2 orp~^. There are also less 
coherent swathes of power at about 0.8 and 4.4 orp~^. The 4.4-orp~^ swath is of 
some interest because it shows that the pattern period of the inner-disk structure 
is still evolving. For m = 3, the same swathes of power that are in the m = 2 
periodogram are present. However, power for the innermost swath is much more 
well-defined, and the pattern period is closer to 4 orp~^; again, there is an indication 
of structure evolution. In addition, a swath of power at the 1.5-orp~^ corotation is 
present in the m — 3 periodogram, but is largely absent in the m — 2 periodogram. 
By m = 4, most of the structure is relatively muted. As in the M2004 disk, the 
commensurate power for modes m — 1, 2, and 3 are likely indicative of some form of 
mode coupling among low-order m, and ILR-corotation-OLR alignments like those 
present in the M2004 periodograms, can be found for m = 2 and 3. 

Even though there are several persistent spiral structures in the disk, the mass 
fiux profile is confusing. Nonetheless, there are some alignments of interest. The 
principal inflow/outflow boundary is at about 15 AU, which is very close to the 



6.7. COHERENT STRUCTURE 183 

corotation of the 4.4- and 4-orp~^ swathes for m = 2 and 3, respectively. The second 
mass inflow/outflow boundary at about 20 AU is commensurate with the well-defined 
m — 2 Fourier mode (2 orp~^), and the third inflow/outflow boundary at about 27 AU 
is commensurate with the corotation of 1.5-orp~^ swath in the m — 3 periodogram. 

To evaluate further the importance of these m-arm structures for mass transport, 
I decompose the BDNL disk into its respective density Fourier components and re- 
construct the disk with only certain modes. Figure 6.16 shows the results. The m — 1 
component has strong power in the periodogram, but it is a minor component of the 
torque profile until r > 30 AU. Most of the power inside r ~ 30 AU is generated by 
spirals with arms m = 2, 3, and 4, and so mass transport should be modeled fairly 
accurately. Regardless, this does emphasize that routines for freeing the star in the 
SV and in CHYMERA need to be developed. 
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Figure 6.14 M2004 disk periodograms showing m = 1, 2, 3, and 4. The curves indicate 
corotation and the ILR and OLR (only OLR and corotation for m = 1). 
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Figure 6.15 BDNL disk periodograms showing m = 1, 2, 3, and 4. The curves indicate 
corotation and the ILR and OLR (only OLR and corotation for m = 1). 
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Figure 6.16 BDNL disk torque reconstruction. The m = 1 mode becomes important 
in the outer disk, but the m = 2, 3, and 4 modes dominate the torque profile for 
r < 30 AU. 
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6.8 Local or Global Transport? 

The locality of mass and angular momentum transort in Gl-active disks has been 
addressed by multiple authors (e.g., Pringle 1981; Laughhn & Rozyczka 1996; Balbus 
& Papaloizou 1999; Gammie 2001; Lodato & Rice 2004; Mejia 2005): Can mass 
transport in Gl-active disks be modeled by a local a prescription, or do the long- 
range torques and/or wave transport that GIs produce make such an approximation 
misleading? If GIs have a global effect on mass transport in the disk, structure 
caused by spatially and temporally variable accretion rates, such as rings or density 
enhancements, can be missed by assuming that the disk evolves like an a disk. 

The above analyses indicate that low-order spiral waves are the dominate non- 
axisymmetric structures, that the mass flux is correlated with resonances of coher- 
ent, persistent low-order spiral waves, and that most of the time-averaged torque 
is produced by these low-order structures. Even though the effective a during the 
asymptotic phase for each disk is similar to what one would expect based on a disk 
theory, I caution anyone from interpreting these results as an indication that one can 
evolve the disk as an a model. The mass fluxes are highly variable, with multiple 
inflow/outflow regions due to the interplay between the large-scale, low-order struc- 
ture. A local model does not seem hkely to represent the mass fluxes in the M2004 
and BDNL disks well. In particular, bursts of GIs would be incorrectly modeled. 
However, if the tcooi proflle were known a priori, which is an unclear prospect, the 



188 CHAPTER 6. GRAVITATIONAL DISK INSTABILITY 

a model would be successful in predicting gross disk properties, such as the typical 
magnitude of the mass flux during the asymptotic phase. 



6.9 Spectral Energy Distribution 

As emphasized by Nelson et al. (2000), an important test of physical relevance for 
a disk simulation is a comparison between SEDs of observed systems and the derived 
SED from a simulation. Because the M2004 simulation was followed through all four 
phases and because the effective temperatures for both simulations are essentially the 
same, I only construct SEDs for the M2004 disk as if it were at some large distance 
d and viewed face-on. For simplicity, I assume that each (r, 0) column of the disk 
emits radiation according to a blackbody law at its effective temperature over its 
surface area. Because I use a cylindrical grid, the ith area element has a solid angle 
dfli — ridrd(l)/(P. The specific fiux, or fiux density, then can be taUied by 



where is the Planck function and Tj is the effective temperature of the ith cell. 
To avoid using distance, I choose to express the SED in terms of AircfuFi,, which is 
a typical observer's SED. As a basis for comparison, I adopt the approach of Nelson 
et al. (2000) and define a fiducial SED to be one that is derived from a disk with a 




(6.5) 
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temperature law r~^-^, which is the median best fit law to the T Tauri disk sample 
presented in Beckwith et al. (1990). 

Figure 6.17 shows the SEDs derived from M2004 for three different stages in its 
evolution. The short dashed line dehneates the SED for the disk near its brightest 
period during the burst phase (2.4 orp), when the luminosity is nearly 19 Lq as 
integrated between 10^° and 10^^ GHz. The long-dashed line dehneates the SED for 
the disk as it enters the asymptotic phase (10 orp), and the solid black line indicates 
the SED for the disk at the end of the calculation (16 orp). The star, which is 
assumed to have an i? = 2Rq and a Tefr — 4000 K, is included in the SED profile. 
The SED has dips in specific luminosity as it nears the star because it has a 2.3 
AU hole and is missing a contribution from an inner and hotter portion of the disk. 
The red lines indicate fiducial SEDs based on assuming a Tefr ~ r~°'^ temperature 
profile for a 0.0033, 0.01, and 0.033 Lq disk and integrating between and 60 AU. 
The fiducial SEDs also have a shght dip in their profile just before transitioning to 
the stellar portion of the SED due to discretizing the temperature profile into grid 
cells 1/6 AU wide in r. The blue line is a fiducial SED calculated for a luminosity 
of 0.0024 Lq, which is the disk luminosity at 16 orp, in the same way as the other 
fiducials but with a 2.3 AU hole. Even though the actual effective temperature profile 
for M2004 is an exponential, the SED for M2004 is observationally consistent with a 
profile Teff ~ ^--o-e because the inner 20 AU of the disk closely follows a T^s ~ r~°'^^. 

Although GIs appear to lead to an exponential Teff profile for this disk, stellar 
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irradiation (D'Alessio et al. 1998, 2001) probably keeps the outer regions of the disk 
warmer than what is modeled here, except for cases of strong shadowing. The ir- 
radiation routine developed by Mejia (2004) is not employed in these calculations 
because, in her routine, the stellar irradiation is directly incident on the disk. The 
optical photosphere for protoplanetary disks is typically three disk scale heights above 
the midplane (D'Alessio et al. 1998, 2001, 2006). Directly modehng this low-density, 
high-temperature atmosphere, where the stellar irradiation is absorbed and scattered, 
along with the relatively high-density, low-temperature disk is a computational chal- 
lenge. 

As will be discussed in more detail in Chapter 7, the SEDs of the M2004 disk also 
suggest that a disk undergoing a burst due to GIs can become very bright and lead 
to very high accretion rates. The M2004 disk reaches a peak luminosity of 19 Lq at 
about 2.4 ORPs, which makes the disk approximately 10,000 times brighter in the 
burst phase than in the asymptotic phase. Furthermore, the disk has a luminosity of 
only about 0.01 Lq around 2.1 ORPs, so the disk luminosity increases by a factor of 
about 1,000 in 80 yr. The increase in luminosity is sudden and may be much shorter 
than reported here due to coarse temporal resolution of the files required to generate 
SEDs. Bursting disks may be related to the FU Ori phenomenon (see Chapter 1). 
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Figure 6.17 The solid line indicates the SED for the M2004 disk at the end of the 
calculation (16 orp or 4,050 yr). The short-dashed and long-dashed curves delineate 
the SEDs at 2.5 orp and at 10 orp. The red curves represent fiducial SEDs for 0.0033, 
0.01, 0.033 Lq for a disk that extends from 1/6 AU to 60 AU and has a T^g ~ r""^ 
for q = 0.6. The blue dashed curve indicates the SED for a disk that is truncated at 
2.3 AU with a g = 0.6. 
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6.10 Conclusions 

The main results of this comparative study are as follows. The disks evolve to- 
ward a thermally self-regulating, gravitationally unstable state, where work and shock 
heating balance radiative losses. The mass transport is highly variable and dominated 
by global, low-order modes. The gross mass transport is similar in magnitude to that 
expected by a local model. However, evolving the disk as an a model will miss the 
large mass flux variabihty and will neglect the spiral shocks and shock bores in a 
disk, which may be important for generating turbulence (Boley et al. 2005) and for 
chondrule formation (Chapter 7). Moreover, the cooling times are too long to lead to 
fragmentation, and as discussed in §5.4.3, convection eddies are unable to lower the 
cooling times to those that result in disk fragmentation. 

What is the fate of these disks? As discussed in §6.6, the asymptotic phase can 
only be maintained for a fewxlO^ yr. Because the models are incomplete, i.e., I 
exclude the inner disk, irradiation, possible infall, dust settling, grain growth, and 
other mass transport mechanisms, it is difficult to comment on the disk's ultimate 
fate. Regardless, I speculate about the big picture in Chapter 8 based largely on the 
results of this comparative study. 



Chapter 7 

DEAD ZONE MODELS 

In Chapter 6, I discussed two disk simulations, where one was evolved with the 
M2004 radiation scheme and the other with the BDNL scheme (Chapter 2). The 
results indicate that (1) disk evolution is sensitive to the details of the treatment 
of radiation physics, (2) cooling times are too long for disk fragmentation, (3) the 
effective a ~ 10~^ during the asymptotic phase, (4) high mass fluxes occur during GI 
bursts (~ 10~^Mq yr~^), and (5) the shocks during the burst rapidly heat the disk 
such that the disk can outshine the primary. The high mass fluxes during bursts may 
indicate an internal, repeatable process for driving the FU Ori phenomenon (Armitage 
et al. 2001). As discussed in Chapter 1, an FU Ori outburst is characterized by a 
rapid (1-lOs yr) increase in optical brightness of a young T Tauri object, typically 5 
magnitudes, with mass accretion rates of the order 10~^Mq yr~^ from the disk onto 
the star (Hartmann & Kenyon 1996). The most promising explanation for the FU 
Ori phenomenon is a thermal instability (Bell & Lin 1994). But even with this model, 
the driving mechanism that initiates the thermal instability remains unclear. In this 
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Chapter, I consider the possibihty that a burst of GI activity in a mass-enhanced 
region of a disk could produce a cascade of instabihties that ultimately leads to an 
FU Ori event. I envision the mass enhancement to be a result of reduced accretion 
in a dead zone, where the MRI-active layer is only a small fraction of the total mass 
in a given column. 

If episodic GI bursts drive the FU Ori phenomenon, then might these bursts 
also produce shocks that are strong enough to form chondrules? lida et al. (2001), 
Desch & Connolly (2002), Ciesla & Hood (2002), and Miura & Nakamoto (2006) 
demonstrated that chondrules can be formed by gas drag friction when a chondritic 
precursor enters a shock. Certain post-shock conditions allow for melts to cool at 
rates consistent with laboratory experiments. One should note that the shock model 
for chondrule formation is not a model for the mechanism driving the shock, and any 
strong shock has the potential to form chondrules. Global spiral waves were identified 
as a potential shock-driving mechanism by Wood (1996). However, it remains unclear 
how these waves can be repeatedly generated. 

For this study, I adopt the hypothesis that bursts of GI activity in dead zones 
drive the FU Ori phenomenon and produce shocks with chondrule-forming pre-shock 
conditions. In order to investigate this scenario, I designed a numerical experiment 
to evolve a massive, highly unstable disk with an initial radial extent between 2 and 
10 AU. The experiment is conducted with three principal objectives in mind: (1) To 
characterize shocks in a disk with a very strong burst of GI activity near 5 AU and 
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to determine whether such shocks can drive chondrule formation, (2) to determine 
whether a strong burst can drive accretion rates high enough to start a cascade of 
instabihties, resulting in an FU Ori outburst, and (3) to investigate coohng times in 
the 2-10 AU region for a massive disk. 

7.1 Expected Shock Strengths 

The shock speed ui for a fluid element entering a logarithmic spiral wave with 
pitch angle i is described by 



where I have assumed that the gas azimuthal motion is Keplerian and where Tp is the 
pattern radius for some global spiral wave. Neglect Vr for simplicity. This leads to 



Notice the limiting behavior of equation (7.2). When Vp 3> r, the shock speed limits 
to the Keplerian speed times sini, and so whether spiral waves with large rp produce 
chondrules is mainly dependent on the pitch angle of the spiral wave. In contrast, 
when r ^ Vp, the speed increases as r, and even shallow pitch angles can produce 
chondrules, depending on the model. 




(7.1) 




(7.2) 
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r(AU) 

Figure 7.1 Expected shock speeds u based on equation (7.2). In the legend, u{l, 10) = 
Uilvp = 1 AU, i = 10°) for pattern radius Vp and pitch angle i. The colored regions 
highlight where chondrule formation is expected based on the Desch & Connolly 
(2002) shock calculations (see text), with the yellow region appropriate for a MMSN 
density distribution and the orange for the same density distribution but with 10 x 
the mass. Shocks that occur inside corotaton for i = 10° will not produce chondrules 
between 1 and 5 AU. However, chondrules can be produced by these low pitch angle 
spirals in shocks outside corotation. If the pitch angle is fairly open, such as i ~ 30°, 
then a spiral wave with a corotation near 5 AU can produce chondrules in the asteroid 
belt and at comet distances. 
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For example, consider the Minimum Mass Solar Nebula (MMSN), where the mid- 
plane p{r) = 1.4 X 10~^ (r/1 AU)"^^^'' g cm~^ (Hayashi et al. 1985). Figure 7.1 shows 
the u-p plane with heavy, solid curves indicating shock speeds for = 1, 3, 5, and 
10^ AU and with i — 10 and 30°. The colored regions on the plot highlight where 
chondrule formation is expected in a MMSN (yellow) and a lOxMMSM (orange). 
The chondrule- forming curves are based on the results of Desch & Connolly (2002), 
with ui ~ — (11 -|- 2 log p (g cm~^)) ±0.5 km s~^. The boundaries set by these curves 
are meant to be illustrative, not definitive. 

Spiral shocks along shallow pitch angle spiral waves (i 10°) are mostly if not 
entirely out of the chondrule-forming range between r = 1 and 5 AU for shocks inside 
corotation. However, spirals with corotation in the inner disk can produce chondrule- 
forming shocks for a wide range of radii, depending on the actual r^. If the spiral 
waves have very open pitch angles (i > 30°), then spiral waves with corotation near 
r — 5 AU can effect chondrule-forming shocks near the asteroid belt and at cometary 
distances. The mass of the disk does not change these general behaviors. 

A major caveat for this simple-minded approach is that the fluid elements' motions 
may be poorly described by equation (7.2) if vertical and radial excursions induced 
by shock bores cannot be neglected. As demonstrated in Chapter 5, strong spiral 
shocks can significantly alter the fiuid fiow. Regardless, this analysis serves as a base 
description for expected shock speeds in the disk for a given spiral wave. Because I 
am interested in exploring the production of chondrules near the asteroid belt and the 
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annealing of solids in the outer disk, I have tried to design a numerical experiment 
that is biased toward producing strong shocks with a corotation near 5 AU and with 
open pitch angles. 



7.2 Methodology 

The numerical experiment is set up as follows. The initial model is the flat-Q disk 
with the surface density enhanced ring (§4.2). The disk is initially between about 

2 and 10 AU in radius, and the total disk mass is 0.167 Mq. The surface density 
ring enhancement is distributed as a Gaussian centered on 5 AU with a FWHM of 

3 AU. The FWHM was chosen to be comparable to the most unstable wavelength 
at 5 AU (§4.2). The mass- weighted Q and average surface density profiles for the 
disk are shown in Figure 7.2. For comparison, the initial E ~ is also shown. 
The instabilties activate near 4.5 AU once the region between 3 and 6 AU reaches a 
mass- weighted average Q ~ 1, even though the minimum Q at 4.5 AU is well below 
unity. 

For the treatment of H2, I consider the arguments in §2.3.3. Recall that energetic 
particles (EP) are attenuated with a surface density of 100 g cm~^. Figure 7.3 com- 
pares the height above the midplane where ■m{z) = pdz — 100 g cm"^ = rric with 
the scale height of the disk, h — m{0) / p. In addition, it also shows the fraction of the 
half disk column mass that is contained above one vertical attenuation length. One 
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should note that EPs will hkely be on random or obhque trajectories, and so EPs may 
affect less of the disk than assumed here. Cosmic ray ionization will be negligible for 
most of the disk inside about 8 AU, and the disk is too cold for thermal ionization 
of alkalis. Production of H+ and will be limited to radioactive decay, with ^^Al 
dominating the EP production at a rate of C ~ 10^^*^ s^^ (Stepinski 1992; Chapter 
2). This low rate will result in a frozen ortho/para ratio for most of the disk. How- 
ever, because there will be some protonated species in the disk, the ratio can evolve 
slowly to its thermal equilibrium ratio while acting out of thermal equilibrium on the 
dynamic timescale. Finally, even for the areas where cosmic rays can penetrate, the 
depletion timescale due to dust is likely shorter than the dynamic timescale in this 
disk everywhere. I use a frozen ortho/para ratio of 3:1 in this study for these reasons. 

The disk is evolved with the BDNL radiation physics algorithm. The grain size 
distribution is assumed to follow the ISM power law, dn ~ a'^-^da (D'Alessio et 
al. 2001), where the maximum grain size Omax is chosen to be 1 mm, and a weighted 
combination between Rosseland and Planck mean opacities is used (see Chapter 2). 
The particles are assumed to be well-mixed with the gas. However, the effects of dust 
settling are heuristically considered by evolving the disk in several ways: standard 
opacity and 1/100, 1/1,000, and 1/10,000 of the standard opacity (Table 7.1). This 
choice in opacities varies the midplane Rosseland mean optical depths from about 
T = 10"^ to unity. It should be noted that the standard opacity simulation has 
Tv/c~ 1 near the inner edge, and so it is close to the limit {tv/c > 1) where advection 
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of photons becomes important (Krumholz et al. 2006), but this effect is not modeled 
in these simulations. For a base simulation, the disk is also evolved adiabatically for 
about 1 orp, where 1 orp for these simulations is one outer rotation period at 10 AU, 
or about 32 yr. In addition, no external radiation is assumed to be shining onto the 
disk, except for a 3 K background temperature. Although such a naked disk is likely 
to be unrealistic, except for a case of strong self-shadowing, the assumption is made 
to bias the disk toward strong shocks inasmuch as irradiation weakens GI strength 
(Cai et al. 2007). If chondrule-producing shocks cannot be created in simulations 
biased in their favor, then it would present a serious problem for GIs as the source of 
chondrule processing. 

Due to inefficient cooling in the standard simulation, high enough temperatures in 
the disk are reached after 80 yr to create radiative timescales that are too short to re- 
solve; the simulation is stopped. As described below, the reduced opacity simulations 
are able to cool much more efficiently, and the disk temperatures remain manageable 
with the time-explicit algorithms used in CHYMERA. However, the reduced opac- 
ity simulations are stopped after about 110 yr due to the development of a strong, 
one-arm spiral, which is treated incorrectly with the fixed star assumption. 
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Figure 7.2 Top: Surface density profile for tlie flat-Q disk with (solid curve) and 
without (dashed curve) the surface density enhanced ring. The surface density has a 
maximum at about 4.5 AU instead of at the target 5 AU due to the steep initial surface 
density profile. Bottom: The initial Q profile with (solid curve) and without (dashed 
curve) the surface density enhancement. Even though the minimum Q drops below 
unity, the GI burst only activates when the mass-weighted average Q approaches 
unity over a 3-AU annulus centered on r = 4.5 AU. 
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Figure 7.3 Left: The scale height of the disk (sohd curve) is shown with the height 
above the disk where EPs are expected to pass through one attenuation length (dotted 
curve). Right: The fraction of the half surface density that is contained within the 
first attenuation length. 
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Table 7.1 Simulation information. The 1024 1/1,000 and 1/10,000 k simulations 
are still being evolved at the time of this writing, and all 1024 simulations will be 
extended to 110 yr if possible for future work. The name of the simulation indicates 
the fraction of the standard opacity that is used during the evolution of the disk. 



Sim. Name 


Resolution r, (j), z 


Duration 


Adiabatic 


256, 512, 64 


33 yr 


512 Standard 


256, 512, 64 


80 yr 


512 1/100 K 


256, 512, 64 


110 yr 


512 1/1,000 K 


256, 512, 64 


110 yr 


512 1/10,000 K 


256, 512, 64 


110 yr 


1024 Standard 


512, 1024, 128 


55 yr 


1024 1/100 K 


512, 1024, 128 


55 yr 


1024 1/1,000 K 


512, 1024, 128 


33 yr 


1024 1/10,000 K 


512, 1024, 128 


33 yr 
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Each simulation is evolved, at least through the first 30 yr, at two resolutions: 
(r, (j), z) = (256,512,64) and (r, 0, z) = (512,1024,128). The lower resolution 
simulation (512 sim) has a grid spacing of 0.05 AU per cell in r and z and the higher 
0.025 AU (1024 sim). For both resolutions, rA0 fa Ar at about r = 4 AU. As 
described below, shock heating is more intense in the 1024 simulations, and so the 
1024 standard is stopped after only about 55 yr for the same reason that the 512 
standard is stopped. The 1/1,000 and 1/10,000 k simulations are also being run at 
high resolution and have been evolved for about 33 yr. For each simulation, 10^ fluid 
elements are randomly distributed in r between about 3 and 7 AU, in over 27r, and 
in z roughly within the scale height of the disk. The fluid elements are integrated as 
the simulation is evolved. 

7.3 Evolution 

Surface density plots for the 512 simulations are shown in Figures 7.4 and 7.5 at 
t ~ 33 and 77 yr. The images at t ~ 33 yr indicate that the burst has a well-defined, 
three-arm spiral and that the spiral arms become denser as the opacity is lowered. By 
77 yr, the disk structures are noticeably different. Each disk has a visually distinct 
two-arm spiral except for the 1/10,000 k simulation. The lower opacity simulations 
appear to have stronger amplitudes, with denser spiral arms and stronger mid-order 
Fourier components (e.g., m = 4 through 10 in Figure 7.7). As will be discussed in 
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more detail below, the 1/10,000 k, simulation is on the verge of fragmentation. 

For comparison with Figures 7.4 and 7.5, the 1024 simulations are shown in Figure 

7.6 at i ~ 33 yr, which is the time when significant deviations in spiral morphology 
become clearly discernible. The higher resolution simulations have additional fine 
structure, but they do not develop clumps or dense knots. In fact, some of the knotty 
structure that develops in the 512 1/10,000 k, simulation is absent at higher resolution. 

As discussed in Chapter 6, visual differences can be quantified by a Fourier mode 
analysis. The {A^) for each disk, averaged over 55-77 yr, is 1.19, 1.43, 2.30, and 2.40 
for the 512 standard, 1/100, 1/1,000, and 1/10,000 k, simulations, respectively. Figure 

7.7 demonstrates the corresponding spectra. The bars denote typical fluctuations, 
and they should not be mistaken for error bars. The low-order structure clearly 
dominates, with the high-order structure falling off steeply. For each simulation, the 
m — 2 mode dominates during this time interval, but m — 1, 3, and 4 are strong 
in all simulations. The 512 standard and 1/100 k, data track each other, with some 
variation in the power in the low- to mid-order Fourier compnents. The 512 1/1,000 
and 1/10,000 k, data also track each other, but diverge from the 512 standard and 
1/100 K amplitudes. In particular, the power in the Fourier modes falls off more 
gently than it does in the 512 standard and 1/100 k, simulations, which as discussed 
above, is noticeable in the surface density plots. 
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Figure 7.4 Surface density maps for the 512 standard and 1/100 k, simulations during 
the burst and shortly after. The strength of the spiral waves increases as the opacity 
is lowered. 
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Figure 7.5 Surface density maps for the 512 1/1,000 and 1/10,000 k, simulations, which 
complement Figure 7.4. 
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Figure 7.6 Surface density maps for the 1024 simulations at approximately the same 
time as in the top panels in Figures 7.4 and 7.5. 
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The energetics of the disk simulations are shown in Figures 7.8 and 7.9, where cu- 
mulative energy loss by radiation (blue) and cumulative energy dissipation by shocks 
(red) are plotted. For the 512 standard case, the shock heating clearly dominates 
over the cooling; the disk is heating up over the entire evolution. Energy is trans- 
ported inefficiently for the highly optically thick disk, and the disk evolves much like 
an adiabatic simulation. When the opacity is reduced by a factor of 100, the cool- 
ing becomes much more important than in the standard simulation, with radiation 
energy losses becoming more similar to heating by shocks. In addition, the shock 
heating rate becomes larger during the burst. When the opacity is reduced by a 
factor of 10^, the radiative cooling is even more efficient and surpasses shock heating. 
The 10^ opacity reduction shows the strongest shock heating and the fastest disk 
cooling. As one would expect from analytic arguments, the opacity has a profound 
effect on disk coohng and, consequently, on spiral shocks. For the 1024 standard and 
1/100 K, simulations, coohng and shock heating are similar to the 512 simulations for 
about the ffist 32 yr, but the 1024 runs exhibit more shock heating, while the cooling 
remains roughly the same. This additional shock heating pushes the disks toward 
GI stability faster than what is seen in the 512 simulations. The departure of the 
1024 shock heating from the 512 curves is commensurate with the visual deviations 
of spiral structure in the surface density images. 
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<Aj^> lor 55-77 yr 




Mode 

Figure 7.7 A„i spectra for the standard, 1/100, 1/1.000, and 1/10,000 k simulations. 
The bars indicate typical fluctuations during the 55-77 yr time period, and should 
not be mistaken for error bars. 
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Figure 7.8 Radiation energy losses and shock heating for the 512 and 1024 standard 
and 1/100 k simulations. The effect of opacity on the energy loss and shock heating is 
clear: the lower the opacity, the faster the cooling. In addition, the 1024 simulations 
exhibit additional shock heating, but comparable cooling to the 512 resolution sim- 
ulations. The short, adiabatic curve is difficult to see because it tracks the standard 
curves closely. 



212 



CHAPTER 7. DEAD ZONE MODELS 




Figure 7.9 Complement to Figure 7.8: Radiation energy losses and shock heating for 
the 512 1/1,000 and 1/10,000 k simulations are shown. 
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The effect of opacity on energy losses in the 512 disks is also demonstrated in 
Figures 7.10 and 7.11. In these figures, brightness temperature maps are shown 
for the same time as the top rows in Figures 7.4 and 7.5, where Tj, — (Trl^/a)^^'^ 
and /+ is vertical outward intensity (see Chapter 2). In addition, I define a cooling 
time temperature Tc = {J^edz/a /o°° — V • FdzY^^, which is the temperature that 
corresponds to a given column's effective fiux if all of the energy were to leave the 
column vertically. If the column is being heated, Tc is set to zero. 

As the opacity is lowered, the spiral structure becomes more clearly outlined, and 
the disk becomes brighter because the photons can leave from hotter regions. The 
brightness maps also demonstrate that sustained fast cooling by convection is absent. 
If hot gas near an optically thick midplane is quickly transported to altitudes where 
T ~ 1, convection can, in principle, enhance disk cooling. However, the efficacy of 
convective cooling is controlled by the rate that energy can be radiated from the pho- 
tosphere of the disk. If localized convective ffows were responsible for fast coohng in 
the optically thick disks, one would expect to find strong, localized T;, enhancements. 
For example, the 512 standard and 1/100 k, simulations would need to have regions 
as bright as the 512 1/10,000 k simulation. The corresponding energy fiux that is 
needed for convection to be sustained and to cool the disk does not occur (see Chapter 
5 for additional discussion regarding convection). Finally, the Tc maps evince that 
radiation can lead to net heating in regions near strong shocks, e.g., the black outlines 
(areas experiencing net radiative heating) around the thin, red spiral waves. 
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Figure 7.10 Tf, and maps for the standard and 1/100 k, simulations. As the opacity 
is lowered, the disks become much more efficient at cooling. See also Figure 7.11. 
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Figure 7.11 Complement to Figure 7.10: Brightness and Tc maps for the 1/1,000 and 
1/10,000 K simulations. 
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The azimuthally averaged brightness temperature profiles are shown in Figure 
7.12. As shown in Figures 7.10 and 7.11, the brightness temperatures are higher 
for the lower opacity simulations. In addition, each profile can be roughly fit by a 
power law ~ r'P, with p = 0.61, 0.84, 0.91, and 0.85 for the 512 standard, 1/100, 
1/1,000, and 1/10,000 k, simulations, respectively. These power laws are consistent 
with observed effective temperature profiles for irradiated disks (e.g, Beckwith et 
al. 1990; Kitamura et al. 2002), but these disks have not reached an asymptotic 
phase, during which the effective profiles may be closer to those found in the M2004 
and BDNL simulations (Chapter 6). 

The importance of the opacity for disk cooling is also shown in Figures 7.13-7.16. 
On these plots, three quantites are shown: the Toomre Q, the mass- weighted Fi, and 
the icooi^//(ri) profiles. For the icooi curves, the cooling times are calculated by 
dividing the azimuthally and vertically integrated internal energy by the azimuthally 
and vertically integrated radiative cooling for each annulus in the disk. The /(Fi) is 
the critical value of tcooi^, below which fragmentation is expected, for the correspond- 
ing Fi. Rice et al. (2005) demonstrated that /(7/5) ^ 12 and /(5/3) 6. Based on 
these values, I assume /(Fi) fti — 23Fi -|- 44 for the stability analysis presented here. 
One should keep in mind that /(Fi) is not a strict threshold and that the relation is 
approximate, especially for simulations that permit the cooling time to evolve with 
the disk (Johnson & Gammie 2003). Regardless, it serves as a general indicator for 
disk fragmentaiton. 
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Figure 7.12 Brightness temperature profiles for tlie standard, 1/100, 1/1,000, and 
1/10,000 K simulations at about 77 yr. Each profile can be fit approximately by a 
power law Ti, ~ r^^ (smooth curves), with p = 0.61, 0.84, 0.91, and 0.85 for the 
standard, 1/100, 1/1,000, and 1/10,000 k simulations, respectively. 
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Figure 7.13 Stability assessment for the 512 standard simulation. The ordinate shows 
values for Q, Ti, and tcooi^/ f(Xi) as a function of r, where /(Fi) ^ — 23ri + 44. If 
tcooi^/ f(Xi) > 1, then the disk is stable against fragmentation. 
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Figure 7.16 Same as Figure 7.13, but for the 512 1/10,000 k simulation. 
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Figure 7.17 Surface density profiles for the 512 standard, 1/100, 1/1,000, and 1/10,000 
K simulations at about 77 yr. Each profile, except for the 1/10,000 k simulation, is 
consistent with a Gaussian between about 2 and 8 AU. The initial surface density 
profile with the density enhancement is shown with +s. The burst reorders each disk 
efficiently. 
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Figure 7.18 Mass fluxes for tlie 512 standard, 1/100, 1/1,000, and 1/10,000 k sim- 
ulations during the 0-55 yr, 33-55 yr, and 55-77 yr periods. For each simulation, 
inward and outward mass fluxes are well above 10"*^ Mq yr~^. Positive mass fluxes 
are inward. 
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As Figures 7.13-7.16 demonstrate, the disk is approaching a state of constant Q 
for a wide range of radii, as expected in an asymptotic phase (Chapter 6). The 
mass-weighted average Qs between 3 and 6 AU are 1.73, 1.69, 1.49, and 1.43 for the 
standard, 1/100, 1/1,000, and 1/10,000 k, simulations, respectively. The variation in 
the average Q is consistent with the measurements. 

Because the cooling times fluctuate rapidly, I average the integrated coohng rates 
and the internal energies for the 65- and 77-year snapshots. The cooling time profiles 
for the appropriate Fi only drop substantially below unity in regions of high Q; these 
disks are stable against fragmentation. As the opacity is lowered, the cooling times 
decrease as well, with the 1/10,000 k, simulation close to the fragmentation limit. It 
should also be noted that the only disk that behaves similarly to a constant Fi disk 
is the standard opacity simulation, and for all other simulations, a constant Fi is 
inappropriate. 

In addition to the changes in the Q profile, the surface density profiles are al- 
tered during the burst, and roughly follow a Gaussian profile between about 2 and 
8 AU (Fig. 7.17), with the 1/10,000 k, simulation showing the most deviation. This 
relaxation to a Gaussian was noted by Boley et al. (2006, 2007c) for the M2004 and 
BDNL simulations, and such a profile seems to be typical of Gl-driven accretion when 
low-order modes dominate. This rapid reorganization of mass also indicates that the 
accretion rates are very high. Figure 7.18 shows mass fluxes in the disk for several 
time intervals. As described in Chapter 6, mass fluxes are calculated by differencing 
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the mass inside a cylinder at two different times. The inward and outward accretion 
rates vary, and can be well above 10""^ Mq yr~^. 



7.4 Shock Strength 

Using the periodogram analysis (Chapter 6), I find that the location of corotation 
for the 171 — 2 and m — 3 spiral waves is at r ~ 4 AU, as measured between 
and 77 yr. The exception is the 512 1/10,000 k, disk, in which weak signals indicate 
corotation near 4 AU for m = 3 and 5 AU for m = 2. In addition, the mass flux 
inflow/outflow boundary in each disk for the 0-55 yr time period (Fig. 7.18) indicates 
that corotation is initially near 4.5 AU. Overall, corotation is about 1 AU inward from 
the target location, but still places chondrule-forming shocks within the asteroid belt 
and at cometary distances for large pitch angle spirals. In all simulations except 
for the 512 1/10,000 k run, the bursts produced spirals with pitch angles of only 
i ft! 10°. Because the 512 1/10,000 k simulation is close to fragmentation (see above), 
the amplitudes for the mid-order Fourier components are large, and the global spirals 
develop a flocculent morphology in areas. Based on the arguments in §7.1, one should 
not expect for many chondrule-producing shocks to be present in these simulations 
unless the vertical and radial motions are important. 
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Figure 7.19 Thermal and spatial histories for a fluid element in the 512 1/10,000 k, 
simulation. Radial and vertical motions vary considerably, as one would expect from 
shock bores. The fluid element experiences many shocks with multiple pre-shock 
conditions. This fluid element's history includes one potential chondrule-forming 
shock at i 20 yr. 
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As discussed in §7.2, fluid elements are used to investigate shock strengths and 
fluctuations in thermodynamic quantities, e.g., temperature, in the 512 simulations. 
Figure 7.19 shows thermal and trajectory histories for a sample fluid element in the 
512 1/10,000 K, simulation that experiences a potential chondrule- forming shock (see 
below). Multiple shocks with a wide range of pre-shock conditions are encountered. 
However, one should keep in mind that shocks in these simulations are unlike sim- 
ple, ID shocks, where the length of the full shock occurs in roughly lO^'^ cm; these 
simulations can only resolve scales of about 10^^ cm or larger. Despite this potential 
difficulty, I use the pressure difference before and after the shock to calculate the Mach 
number Ai. Once A4 is known, the velocity normal to the shock ui, in the frame 
of the shock, can be inferred. Using the pressure histories for each fluid element, a 
possible shock is identified whenever the dp/dt changes from negative to positive then 
back to negative. The pre-shock flow is then taken to be the first sign switch, and the 
post-shock flow is the second sign switch. If Ai'^ — (7 -|- 1) 77/ (27) -|- 1 > 2, the event 
is counted as a shock, where rj = {p2 — P2) /pi is the fractional pressure change and 
7 is the average of the pre- and post-shock first adiabatic index. Once A4 is deter- 
mined, Ui — Aicgi is calculated, where c^i is the pre-shock adiabatic sound speed of 
the gas. Other methods were tried, e.g., estimating the speed of the spiral wave, but 
local variations in pattern speed and three-dimensional shock orientations render such 
methods exceedingly difficult. As can be seen in the pressure panel in Figure 7.19, 
using 77 does a decent job of identifying shocks. Some counting errors and poor esti- 
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mations of shock strength are associated with this technique. In particular, because 
shocks in these simulations are spread out over multiple cells by the AV algorithm, 
post-shock pressures could be reduced by radiation transport and three-dimensional 
wave effects, e.g., shock bores, which would cause underestimates for A4. 

Table 7.2 indicates shock information as extracted from the fluid element informa- 
tion. Column two shows the total number of detected shocks (TS) with Al^ > 2, and 
column three displays the average number of shocks per fluid element. Columns four 
(TS Mass) and five (CS Mass) show the total dust mass (see below) that goes through 
a shock with M"^ > 2 and the total dust mass that encounters a chondrule-forming 
shock, respectively. Finally, column six shows the total number (Total FE) of fluid 
elements that remain after the 70 yr time period, i.e., the elements that were not 
accreted onto the star or fluxed into background density regions. 

As described in §7.1, I consider a shock to be chondrule-forming when ui lies in a 
1 km s~^ band between 5 km s~^ and 11 km s"^ and the pre-shock density is between 
logp( g cm~^) = -8.5 and -10.5. Chondrule-forming shocks should be thought of as 
having the potential to form chondrules, but may not yield chondritic material due 
to incompatible dust to gas ratios, too high or low cooling rates, and fractionation. 
Likewise, by dust processing, I only mean that dust may be altered by the shock. 
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Table 7.2 Shock information for the time period between 10 and 70 yr. TS indicates 
the total shocks encountered by all fluid elements. TS/FE gives the average number of 
shocks for a fluid element. TS Mass is a rough estimate of the total dust mass pushed 
through shocks in each disk, and CS Mass is the total dust mass that encounters a 
chondrule-forming shock. Finally, Total FE indicates the number of fluid elements 
that remain in the simulated disk at the end of the 60 yr time period. The errors, 
when specified, are taken to be Poissonian, and they are only given for comparing the 
simulations to each other. 
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For the processed dust estimates, I assign each fluid element a mass by calculating 
the total disk mass within some Ar and distributing that material evenly among all 
fluid elements in that Ar. In addition, I assume a gas to sohds ratio of 100 everywhere 
for each simulation. I note that this is inconsistent with the assumption that the 
opacity is lowered in three of the simulations because of settling, which would make 
the dust to gas ratio much higher in the midplane and lower everywhere else. Because 
only an order of magnitude estimate is sought, this detail is ignored inasmuch as the 
midplane will process more solids per shock and the high-altitude shocks will process 
less. 

I check whether the total number of detected shocks is reasonable by evaluating 
the expected number of shocks in a Keplerian disk for m-arm spiral waves. If the 
pattern radius for the spirals is rp, the number of shocks during some time period At 
for a fluid element orbiting at r is 



Evaluating equation (7.3) for 1000 fluid elements evenly distributed in annuli between 
2 and 8 AU, with rp — 4 AU, with m — 3, and with At — 60 yr yields approximately 
8000 shocks. So the number of detected shocks in these simulations is reasonable. 

To determine whether the differences in the number of shocks between each simu- 
lation are significant, I take the uncertainty in these measurements to be Poissonian. 




(7.3) 
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Within the uncertainty, the number of shocks per fluid element is roughly consistent 
for each disk except for the 512 1/10,000 k run. The more flocculent spiral structure 
in this simulation not only produces a larger number of shocks in the disk, but also 
effects chondrule-producing shocks. Several thousand of dust are pushed through 
shocks in each simulation, with a few M® of dust experiencing chondrule-forming 
shocks in the 512 1/10,000 k, run. I remind the reader that these are only order of 
magnitude estimates. 

All chondrule-forming shocks transpire within the first 40 yr of evolution, with 
three occurring within the first 21 yr (Table 7.3), and so the onset of the burst creates 
the strongest spiral shocks. These shocks all occur between 3 and 5 AU. Because the 
pattern radii for the global spiral waves are between 4 and 5 AU, radial and verti- 
cal motions and/or unusual wave geometry are likely responsible for the high shock 
speeds. Figure 7.20 shows shock events with M."^ > 2 on the Ui-p plane. The green 
strip indicates the chondrule-forming region. Notice that there are a number of other 
shocks that fall just outside the green area. The strip is only a rough estimate based 
primarily on one group's ID shock simulations that are known to require refinements. 
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Table 7.3 Information for chondrule-forming shocks in the 512 1/10,000 k simulation. 
Columns three and four indicate the r and z where the shock occurs, while columns 
five and six indicate the initial r and z for each fluid element. All chondrule-forming 
shocks take place within the first three vertical cells {Az — 0.05 AU), which is also 
within roughly a third of the local disk scale height for these r. 
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In addition to chondrule formation and dust processing, the relative importance 
of a given shock strength in heating the disk can be evaluated. Figure 7.21 shows the 
number of shocks N per AM. on a log-log plot. Although there is variation among 
the profiles, each simulation roughly shows NAM. ~ M.''^AA4 for p f» 4. With this 
in mind, consider a specific shock energy e ~ A4^. The energy added to the disk per 
AM is then EAM ~ M^NAM. Based on this argument, dissipation is only evenly 
distributed among the shocks for NAM ~ M~^AM. Even though the shocks are 
mediated by global spiral waves during these bursts, the greatest contribution to the 
heating by shocks comes from weak shocks over the r-range for which fluid trajectories 
have been followed. 

7.5 Discussion and Conclusions 

In this section, I review the implications of the results of this study for disk frag- 
mentation, the effects of opacity on disk cooling, the FU Ori phenomenon, and chon- 
drule formation. 1 remind the reader that the simulations presented here are meant to 
be a numerical experiment that explores the possible connection between chondrules, 
FU Ori outbursts, and bursts of Gl-activity. Moreover, this experiment provides a 
systematic study of the effects of opacity on disk cooling, which complements the Cai 
et al. (2006) metallicity study and provides a test bed for disk fragmentation criteria. 
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{km s" ) 

Figure 7.20 Shocks on the Ui-p plane for the 512 1/10,000 k simulations. The green 
area shows the chondrule-forming region, based mainly on the results of Desch & 
Connolly (2002) (see §7.1). The green strip should not be thought of as definitive, and 
slight changes in the strip's location could identify more chondrule-forming shocks, 
particularly at low ui. 
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Figure 7.21 The number of shocks versus Ai, binned in AAi 
profiles roughly follow a power law. 
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7.5.1 Fragmentation 

After the onset of the burst, none of the disks fragments, and only the 512 1/10,000 
K simulation shows dense knot formation during the peak of the burst (10-30 yr). 
These knots do not break from the spiral wave even in the 1024 1/10,000 k simulation, 
and so clump formation does not seem to be missed due to poor resolution. One should 
also keep in mind that for this simulation, the opacity is suddenly dropped by a factor 
of 10^, and such knot formation may not occur in a disk with more realistic settling 
timescales. As discussed below, this is also a caveat for the chondrule formation 
results. On the other hand, it does indicate that disk fragmentation may occur inside 
10 AU under very extreme conditions. 

The stability of these disks against fragmentation is supported by Figures 7.13- 
7.16. These figures indicate that coohng rates for the 512 standard and 1/100 k, 
simulations are too low to cause disk fragmentation. The 512 1/1,000 and 1/10,000 
K, simulations do have areas where ( — tcooi^//(ri) < 1, but the disks approach 
stability against GIs in those regions {Q > 1.7; Durisen et al. 2007a). As noted in 
§7-3, C ^ 1 is not a strict instability criterion, but it does serve as an estimate for 
disk stability against fragmentation. For no simulation is C well below unity where Q 
is also < 1.7. 
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Lowering the opacity increases the coohng rates. The 1/10,000 k, simulation ex- 
hibits the most rapid coohng because the midplane optical depths are near unity, 
which results in the most efficient radiative cooling possible. Changes in dust opacity 
have a profound effect on disk cooling, contrary to what is claimed by Boss (2001, 
2004). Although not modeled, if the opacity were to continue to drop such that the 
midplane optical depth becomes well below unity, cooling would once again become 
inefficient. In addition, supercooling of the high optical depth disks (standard and 
1/100 K.) by convection does not occur. These findings are consistent with analytic 
arguments by Rafikov (2005, 2007) and with numerical studies of disk fragmentation 
criteria by Gammie (2001), Johnson & Gammie (2003), and Rice et al. (2005). 

I have demonstrated that the BDNL radiation algorithm used for this study cou- 
ples the high and low optical depth regimes well (Chapter 3). My conclusions regard- 
ing fragmentation are based on mulitple analyses: surface density rendering (Figures 
7.4-7.6), an energy budget analysis (Figures 7.8 and 7.9), cooling temperature and 
brightness maps (Figures 7.10 and 7.11), and a cooling time stability analysis (Fig- 
ures 7.13-7.16). Furthermore, the resolution in each direction for all simulations was 
doubled, and these 1024 simulations were evolved to at least the peak of the burst 
(t 20-33 yr). Radiative cooling and shock heating curves for the 1024 standard and 
1/100 K simulations were compared with curves for the 512 simulations. The 512 and 
1024 curves show similar radiative losses, but the 1024 simulations exhibit additional 
shock heating. Fragmentation was not missed. 
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It is also pertinent to evince that the code can detect fragmentation when coohng 
rates are high and Q is low. Figure 7.22 shows snapshots for the 512 1/10,000 k 
simulation and a similar simulation but with the divergence of the fluxes artificially 
increased by a factor of two. In the normal simulation, kinks in the spiral waves 
form during the onset of the burst. As discussed above, the disk is very close to the 
fragmentation limit, but the knots do not break from the spiral wave. Figure 7.16 
suggests, although for a later time, that increasing the cooling rates by a factor of 
two should drop ( well below unity in \ow-Q regions. As expected, the 512 1/10,000 
K disk with enhanced cooling fragments. 

Three clumps form, one for each spiral wave, between 4 and 5 AU. The location 
of clump formation is consistent with the prediction by Durisen et al. (2007b) that a 
spiral wave is most susceptible to fragmentation near corotation. One of the clumps 
survives for several orbits and eventually passes through the inner disk boundary. 
Resolution is always a concern for simulations. Because these results are consistent 
with analytic fragmentation limits and numerical fragmentation experiments, it ap- 
pears that CHYMERA can detect fragmentation at the resolutions employed for this 
study. 
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Figure 7.22 The 512 1/10,000 k simulation and a similar calculation, but with the 
cooling rates artificially increased by a factor of 2. Based on Figure 7.16, one would 
expect for the disk to fragment with this cooling enhancement, and it does. 
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7.5.2 FU Ori Outbursts 

In these simulations, the strong bursts of GI activity provide high mass fluxes 
(M > 10~^ Mq yr~^) throughout each disk. Even though corotation is at r ~ 4 AU 
for the major spiral arms, the 2 AU region of the disk is strongly heated. It is not 
difficult to speculate that if a larger extent of the disk were modeled, the temperatures 
due to shocks would be large enough to thermally ionize alkalis and drive a thermal 
MRI. Figure 7.23 indicates maximum, average, and minimum temperatures of the 
fluid elements in the 512 disks. The standard and 1/10,000 k simulations have the 
hottest fluid elements because the gas cannot cool quickly or the shocks are extreme, 
respectively. The other simulations show strong temperature variations as well, but 
the peak temperatures of the disk are not as high. From these simulations it appears 
to be plausible, but by no means proven, that a burst of GI activity as far out as 4 AU 
can activate a thermal MRI inside 1 AU, which may then be responsible for a thermal 
instablility. Although I am simulating a very massive disk, I remind the reader that 
such density enhancements may be plausible for massive T Tauri systems (see Chapter 
1). Additional studies need be conducted, preferably with a self- consistent build-up 
of a dead zone, to address the efficiency of this mechanism in low mass disks. 
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There are at least three observable signatures for this mechanism. First, if a GI- 
bursting mass concentration at a few AU ultimately results in an FU Ori phenomenon, 
then one would expect to see an infrared precursor, with a rise time of approximately 
tens of years. Second, one would also expect for a large abundance of molecular 
species, which would normally be frozen on dust grains, to be present in the gas 
phase during the infrared burst due to shock heating (T. Hartquist 2007, private 
communication). Third, approximately the first ten AU of the disk should have large 
mass flows if the burst takes place near 4-5 AU. I speculate based on these results 
that if the burst were to take place at 1 AU, then high outward mass fluxes should 
be observable out to a few AU. 

7.5.3 Dust Processing 

Each simulation shows that a large fraction of material goes through shocks with 
M"^ > 2. Although these shocks are weak, their abundant numbers may result in 
the processing of dust to some degree everywhere in the 2-10 AU region. In fact, 
such processing may be necessary for prepping chondritic precursors for strong-shock 
survival (Ciesla 2007, private communication). 
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Figure 7.23 Initial r and minimum, average, and maximum temperature plots for 
fluid elements in each 512 case. The temperature variations are quite large, and 
the standard and 1/10,000 approach maximum temperatures of 1000 K. The 

values on the abscissa are due to the fluid elements that were lost from the simulated 
disk. 
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Based on the arguments presented in §7.1, the intent was to produce spirals with 
large pitch angles by constructing a disk biased toward a strong, sudden Gl-activation 
near 5 AU. Even with this bias, the pitch angles of the spirals remain small, with 
i 10°. Why are the pitch angles so small? According to the WKB approximation, 
coti =1 krv/m I (Binney & Tremaine 1987), where kr is the radial wavenumber 
for some m-arm spiral. As discussed in Chapter 4, the most unstable wavelength 
is roughly ~ 27:03/ Qk ~ ^nh/Q for disk scale height h. For a disk unstable 
to nonaxisymmetric modes, corresponds to some m-arm spiral (e.g., Durisen et 
al. 2007b). By relating kr = 27r/3m/A„, 



in the linear WKB limit (| kr-r/m |^ 1), where /? is a factor of order unity. Be- 
cause (3Qr/h ~ 10 in gravitationally unstable disks, the hnear WKB analysis may 
be marginally applicable. Equation (7.4) predicts that linear spiral waves in these 
disks should have pitch angles i fsi 6 to 11° for /3 between 1 and 1/2, respectively. It 
appears that this estimate for i extends accurately to the nonlinear regime. 

These spirals are efficient at heating the disk and transporting angular momen- 
tum, but not at producing chondrules. Only for the 512 1/10,000 n simulation 
are chondrule-producing shocks detected; the 1024 simulations are still being ana- 
lyzed. This low opacity disk is on the verge of fragmentation, and so the presence of 



cot i f3Qr I h 



(7.4) 
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chondrule-forming shocks is likely due to the flocculent spiral morphologies, including 
kinks in spiral waves. 

To estimate the occurrence of a chondrule-forming shock, I employ a generous Ui-p 
criterion. I do not take into account the optical depth criterion of Miura & Nakamoto 
(2006) on grounds that the large scale over which these shocks take place can allow 
for chondrules to equilibrate with their surroundings (Cuzzi & Alexander 2006; see 
Chapter 1). However, one should be aware that the optical depth criterion of Miura 
& Nakamoto will hkely exclude all chondrule-producing shocks inside 4 AU in the 
Ui-p plane for these simulations. All chondrule-forming shocks occur between r — 3 
and 5 AU and at altitudes that are roughly less than a third of the gas scale height 
(Table 7.3). Because a large degree of settling is assumed, these shocks are located 
in regions that may be consistent with the dusty environments in which chondrules 
formed (Wood 1963; Scott & Krot. 2005). 

I estimate that ~ 1 of dust is processed through chondrule-forming shocks. 
Because these shocks are limited to the onset of the Gl burst, a few xlO of dust 
would be processed in the 1/10,000 k, disk if it went through about ten outbursts. If 
more outbursts take place than those that lead to an FU Orionis event, then more 
chondritic material could be produced. One should also remember that to produce 
these shocks, the disk was pushed toward fragmentation by suddenly dropping the 
opacity by a factor of 10,000. In order for bursts near 4 or 5 AU to produce chondrules, 
the disk must be close to fragmentation. 
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An alternative to the fragmentation scenario is suggested in Figure 7.1. In addi- 
tion to driving the FU Ori phenomenon, bursts that occur near 1 AU can produce 
chondrules even with small spiral pitch angles. Moreover, bursts that activate slightly 
further out, near 2-3 AU, could drive chondrule formation beyond the frost line and 
produce the annealed dust that is observed in comets (Harker & Desch 2002; Wooden 
et al. 2005; McKeegan et al. 2006). There are several major advantages to this sce- 
nario, which I discuss in more detail in Chapter 8. 



Chapter 8 



CULMINATING CONCLUSIONS 

In this dissertation, I have addressed the behavior of spiral shocks in protoplane- 
tary disks, and have attempted to make connections between these shocks and mul- 
tiple disk phenomena and planet formation. Spiral waves exhibit rich dynamics, and 
in general, they are by no means simple, two-dimensional beasts. They are likely 
responsible for some of the most energetic accretion events in disks that surround 
newly-forming. Sun-like stars. I find that direct giant planet formation by disk in- 
stability requires demanding and possibly unphysical conditions. Despite this, these 
waves are not inconsequential to building planets. Spiral shocks can serve as the 
machinery for thermally processing dust throughout a planet-forming disk; therefore, 
they play a weighty role in the formation of planets. In this final Chapter, I summa- 
rize the work and the principal conclusions for the studies in this dissertation. I also 
outline the hypothesis that gravitational instabilities in dead zones are responsible 
for chondrule formation and the FU Orionis phenomenon. 
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8.1 Three-dimensionality of Spiral Waves 

As a spiral wave develops into a shock, the post-shock region is typically thrown 
out of vertical force balance and rapidly expands. One should only expect a spiral 
shock to behave like a density wave (two-dimensionally) when the gas is isothermal 
and when self-gravity can be neglected. If self-gravity is important in an isothermal 
shock, then the gas will compress. When the gas is stiffer than isothermal, the gas 
will almost always expand after the shock, except for cases of low Mach number and 
strong self-gravity. This rapid expansion of the gas is called a shock bore, and these 
bores can create breaking waves and induce sudden radial as well as vertical motions 
in the gas. Shock bores may also be an important source of disk turbulence. As 
demonstrated in Chapter 5, large vortical flows are effected by spiral shocks, but 
with the analyses used in these studies, no distinction between large-scale stirring or 
mixing can be made. 

As cooling times decrease, shock heating rates increase and mid-order spiral struc- 
ture becomes stronger. Unless the cooling times drop below the fragmentation limit 
(see below), radiation losses balance shock dissipation and work done by gravity on 
the gas. Shocks are not limited to the midplane, but extend to all altitudes along 
the spiral wave front and elsewhere at mid- to high-disk altitudes as a result of waves 
generated by shock bores. 

The pitch angles of spiral waves in Gl-active disks are relatively small, even for 
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large bursts. Only when the disk approaches fragmentation are pitch angles with 
i > 10° present. Based on fluid element histories, the frequency of shocks as a 
function of Mach number decreases approximately as a power law during bursts of GI 
activity, and the low Mach number shocks contribute the most to heating the disk. 
Shock frequencies during an asymptotic state have yet to be explored. 

As shown in Chapters 6 and 7, low-order structure is primarily responsible for 
mass transport in Gl-active disks, but as the cooling times approach the fragmen- 
tation limit, mid-order structure adds a significant component to the gravitational 
torque. During the asymptotic phase of disk evolution, the effective a approaches the 
predicted a for a local model, despite the mass transport being dominated by global 
modes. Even when the a model is reasonably consistent with GI mass transport 
during the asymptotic phase, the picture of mass slowly diffusing through the disk is 
inaccurate; local mass fluxes in these disks are highly variable. Moreover, using a lo- 
cal model to describe mass transport and disk heating by GIs during the asymptotic 
phase requires knowledge of the cooling rates a priori, which presents a challenge. 
Finally, a local formalism for mass transport fails to capture GI bursts well, even 
though the pitch angles are small, | kr-r/m |~10, for the spirals in these simulations 
(§7.5.3). In Chapter 7, it was shown that bursts of GIs centered on r ~ 4-5 AU cause 
high mass fluxes over the entire simulated disk. This is a concern for implementing a 
local model because bursts are a real possibility in disk evolution, e.g., in producing 
the FU Orionis phenomenon. 
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8.2 Disk Fragmentation 

The radiation hydrodynamics simulations presented in Chapters 6 and 7 have 
been analyzed using a suite of analyses that test multiple aspects of the physics that 
could affect disk fragmentation. In particular, the results address the importance of 
convection, opacity, and radiation boundary conditions (see also Nelson 2006) on disk 
cooling. 

Convection: The convection test results show that both the SV and CHYMERA 
allow for convection under appropriate conditions (§3.3). As discussed in Chapters 5 
and 6, convection is active in the M2004 disk during the axisymmetric phase. While 
active, the convective cooling rate is comparable to the radiative cooling rate, as is 
expected because energy must ultimately be radiated away at the disk's photosphere. 
Once GIs activate, convection is disrupted, and no fragmentation occurs when either 
the M2004, C2006, or BDNL algorithm is used for coohng. 

Mayer et al. (2007) suggest that the M2004 simulation does not portray efficient 
cooling by convection because the optical depths are not large. However, convection- 
like motions are also noted in the flat-Q disk without the density enhancement when 
the opacity law is tuned to induce convection. The cooling times in this massive, 
highly optically thick disk are much longer than those expected to lead to fragmen- 
tation. Moreover, in the dead zone models presented in Chapter 7, brightness tem- 
perature maps were shown to look for local regions where fast cooling by convection 



250 CHAPTER 8. CULMINATING CONCLUSIONS 

might be taking place. No anomalous bright spots were noticeable. In the series of 
dead zone calculations, Rosseland mean midplane optical depths spanned between 
T = 10^ and 1. Supercooling by convection was absent. This is consistent with ana- 
lytic arguments by Rafikov (2007), with our numerical tests, and with work by Nelson 
et al. (2000), who assumed efficient convection in their estimations for the local disk 
effective temperature. Sudden vertical motions in disks along spiral shock fronts are 
shock bores. 

Opacity: A systematic study of the effects of the opacity on the strength of spiral 
waves has been presented in Chapter 7. The study demonstrates that opacity can 
have a profound effect on disk cooling and on the structure of spiral waves. This is 
evinced in the cumulative energy losses, the Fourier structure analysis, local cool- 
ing time plots, and the brightness temperature maps. These results complement the 
study by Cai et al. (2006), but are in disagreement with Boss (2002, 2004a), who 
posits that metallicity does not change the outcome of disk instability due to efficient 
cooling by convection. The studies presented in this dissertation indicate that radia- 
tive cooling ultimately controls convection, which precludes cooling rates that lead to 
disk fragmentation in all disks I have studied (r < 40 AU). I find, as expected from 
analytic arguments, cooling becomes more efficient as the Rosseland mean midplane 
optical depth approaches unity. 

Cooling times: As the opacity is lowered, the cooling times approach the fragmen- 
tation limit. Although the 1/10,000 k, simulation is on the verge of fragmentation 
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and knots do form in the spiral waves during the onset of the burst, none of the spi- 
rals or knots breaks into a clump. When the resolution is doubled in each direction, 
which provides a very high resolution simulation with (r, 0, z) = (512, 1024, 128), the 
knots still do not turn into fragments. The behavior of each simulation is consistent 
with the fragmentation criterion tcooi^//(ri) = C ^ 1) ^-^d the 1/10,000 k, disk only 
fragments when its cooling is artificially enhanced by a factor of two. 

As discussed above, when the opacity in the mass enhanced flat-Q disk is abruptly 
decreased by a factor of 10^, the disk shows fragmentation-like behavior. A similar 
effect is reported by Mayer et al. (2007), who find that their disk only fragments when 
the mean molecular weight is suddenly switched from /j, = 2.4 to = 2.7. It should 
be noted that the M2004 simulation was accidentally run with // = 2.7, and Cai et 
al. (2006) purposefully ran their simulations at the high /i for comparison with the 
M2004 results. Neither of the studies reported disk fragmentation. This may indicate 
that when a disk fragments shortly after a sudden switch in a numerical parameter, 
e.g., opacity here and jj, in Mayer et al., the fragmentation may be numerically driven 
rather than physically. Regardless, the results do indicate that disk fragmentation by 
GIs may be possible under extreme conditions. Demonstrating that such conditions 
can be realistically met remains to be seen. Based on these results, disk fragmentation 
for r < 40 AU appears to be the exception rather than the norm. 

It should also be noted that recent results by Stamatellos et al. (2007), who find 
that disks can fragment at large radii (r > 100 AU), are not contradicted by this 
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study or by analytic arguments (Rafikov 2005). 



8.3 Unified Theory 

A hypothesis behind the work presented in Chapter 7 is that bursts of GI activity, 
dead zones, the FU Ori phenomenon, and chondrule formation are hnked. The general 
picture is that mass builds up in a dead zone due to layered accretion and that GIs 
erupt once Q becomes low. This activation of the instability causes a sudden rise in 
the mass accretion rate that heats up the disk inside 1 AU to temperatures that can 
sustain a thermal MRI. The thermal MRI shortly thereafter activates the thermal 
instability (Armitage et al. 2001; Zhu et al. 2007). In addition, the strong shocks 
process dust and form chondritic material in the asteroid belt and at comet distances. 

The findings presented in Chapter 7 indicate that bursts as far out as about 
4 AU can cause FU Ori mass accretion rates and significantly heat the inner disk. 
Observational signatures of this scenario include high mass fluxes out to about 10 AU, 
an overabundance of typically ice phase molecular species in the gas phase, and a mid 
to far infrared precursor to the optical FU Ori outburst. For all opacities studied, 
multiple fluid elements reached peak temperatures near 800 to 1000 K between about 
2 and 3 AU. The standard opacity simulation (midplane r ~ 10^) and the 1/10,000 
K simulation (r ~ 1) show the highest maximum temperatures. 

The spiral shocks for all simulations maintained relatively small pitch angles {i ~ 
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10) during the burst. The exception is the 1/10,000 k, simulation, which is also the 
only simulation to contain chondrule-forming shocks. This particular run is the closest 
to fragmentation, and has the most flocculent spiral structure. Prom these results, I 
conclude that in order for bursts near 4 or 5 AU to bear chondrule-forming shocks, 
the disk must be on the verge of fragmentation (see also Boss & Durisen 2005a,b). 
Recall that the opacity in 1/10,000 k, disk was suddenly lowered by a factor of 10^ to 
model heuristically an extreme case of rapid growth and settling of solids. As noted 
above, pushing a disk toward fragmentation may only be possible in extreme, and 
possibly unphysical, situations. 

As a result, I return to Figure 7.1. Based on the simple, analytic argument, chon- 
drule formation can take place at cometary distances and in the asteroid belt for 
GI- bursting dead zones between about 1 and 3 AU. As suggested by Wood (2005), 
chondritic parent bodies may be representative of temporal as well as spatial forma- 
tion differences. In particular, he suggests that ordinary chondrites formed near the 
asteroid belt, while carbonaceous chondrules formed beyond the frost line based on 
petrological arguments. Bursts in the inner disk seem to be able to accommodate these 
observations, which is congruent with driving the FU Ori phenomenon (Armitage et 
al. 2001). As the disk evolves, the location of the dead zone can vary, with multiple 
bursts occurring at a few AU. Accretion rates between 10~® and 10~^ Mq yr~^ could 
build up a 0.01 Mq dead zone every lO'' to 10^ years, respectively. Some of these 
bursts may drive the FU Ori phenomenon and produce chondrules, but others may 
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only be responsible for chondrule-formation events. In this scenario, the disk need 
not fragment, and there should be between 10 and 100 separate chondrule-formation 
events. 



8.4 Numerical Improvements 

I have described the most up-to-date implementations of CHYMERA and the 
Standard Version (SV) of the Indiana University Radiation Hydrodynamics Code, 
including various code improvements that were necessary to accomplish the scientific 
objectives of this dissertation work. CHYMERA is different from the SV in one 
significant way: the rotational states of molecular hydrogen are taken into account 
for equilibrium statistics or for a given ortho/para mixture. This modification requires 
that e, not e^/^, be fluxed in the energy equation, which in turn forces the — pV • v to 
be exphcitly calculated during sourcing. The ideal gas equation of state is cast in the 
form p — kpT{e, p) / /imp, where T is calculated from a table based on the cell's specific 
internal energy e = e/p. Dissociation and ionization of gas phase species are ignored 
for the current implementation, but could be added at a later date by changing the 
look-up table to account for a non-constant p. Until this is done, the code can only 
model gas phase hydrogen behavior up to temperatures of approximately 1400 K. 

Explicit modeling of the rotation behavior of molecular hydrogen is necessary 
due to the strong dependence of gas dynamics on the first adiabatic index Fi. For 
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example, Fi changes the behavior of shock bores in disks (Chapter 5) and sets the 
critical cooling time tcooi^ = fi^i) ~ — 23ri + 44 for fragmentation. Moreover, 
seemingly innocent approximations for the first adiabatic index, e.g., Fi = constant 
or quadratic, may cause erroneous dynamics (§§2.3.2 and 7.3). 

In Chapter 2, I reviewed the M2004 and C2006 radiation algorithms, and I in- 
troduced a new routine (the BDNL algorithm) that uses vertical rays to couple the 
optically thin and thick regions of the disk. The accuracy of these algorithms has 
been evaluated (Chapter 3) and shown to follow the expected behavior in a steady 
state test, a contraction test, and a convection test. The limits of each routine have 
also been demonstrated. Each algorithm needs to resolve the disk vertically with 
about five cells when the midplane optical depths are large. Although this number is 
not absolute, it provides a rule of thumb for SV and CHYMERA users. 

The BDNL routine is an improvement over the M2004 and C2006 schemes because 
it allows for complete cell-to-cell coupling in the vertical direction. However, signifi- 
cant improvements can be made by including rays to account for radiation transport 
in all directions. Such an approach would also help to correct the odd temperature 
contours that can arise if diffusion is used in low opacity regimes (§3.4). 
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8.5 Future Code Development 

I present here several suggestions for algorithms that may improve CHYMERA's 
efficiency and extend its repertoire of physics. The suggestions are only meant to 
sketch the possible modifications. These items will be listed on CHYMERA 's Den 
(http://hydro.astro.indiana.edu/chymera), a wiki that provides code updates, docu- 
mentation, and bug information for CHYMERA. Some of the items listed below have 
already been added to the wiki, and some of the suggestions are being pursued by 
myself or other principal CHYMERA users. 

8.5.1 Physics Modifications 

The star should be freed and allowed to respond to disk so that m — 1 spirals and 
subsequent dynamics can be treated well. It is my opinion that this algorithm should 
be implemented before any other improvements are made. 

Using vertical rays to solve for part of the divergence of the flux demonstrated that 
complete cell-to-cell radiative coupling has a noticeable effect on the evolution of the 
disk (Chapter 6). Indirectly, this suggests that such couphng in the other directions 
can have just as important of an effect. The next step toward improving the BDNL 
algorithm could be to include rays along the x and y directions by interpolating 
between cylindrical and Cartesian grids. The z direction can be treated as done in 
the BDNL algorithm. 
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The current opacity calculation in CHYMERA for the BDNL radiation algorithm 
uses the midplane Rosseland optical depth to set the interpolation between the Planck 
and Rosseland means. A way to smoothly interpolate between the Planck and Rosse- 
land mean opacities along the ray is desirable. This will ensure that the Planck mean 
is used in regions of low optical depth regardless of the Rosseland midplane optical 
depth. 

As discussed in §2.4.3, limiters are employed to help amalgamate hydrodynamics, 
shock physics, and radiation physics in the SV and in CHYMERA. The use of limiters 
might be significantly reduced by subcycling the radiation algorithm once the radi- 
ation timescale becomes smaller than, say, 1/100 the Courant time. The radiation 
timescale can be evaluated at the beginning of each step by looking for the minimum 
over the grid of e^/ | V • F |j, which is the ratio of the internal energy density to the 
divergence of the flux for the iih cell. The time step can then be set to some fraction 
of that minimum. 

Stellar irradiation can be employed using the technique described in D'Alessio et 
al. (1999), where stellar irradiation's contribution to the mean intensity is included 
by modeling the optical surface. This method requires knowing the structure of 
the optical photosphere of the disk, which is typically at altitudes three times the 
pressure scale height and which is not modeled in these simulations. However, several 
assumptions can be made to allow one to proceed. For example, a direct correlation 
between the optical surface and the surface at one disk scale height can be assumed. 
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or no correlation between the surfaces can be assumed. Numerical experiments would 
determine the effects of these assumptions, and studies of shock bores in disks with 
extended, hot atmospheres could indicate relations between the disk midplane and 
the optical surface. 

CHYMERA can be modified to evolve a two-fluid grain distribution. As long as 
the dust grain sizes are in the Epstein drag hmit (WeidenschiUing 1977), a distribution 
of grains could be evolved as a pressureless fluid (e.g., Johansen & Klahr 2005). In 
combination with a well-mixed grain component that always remains entrained with 
the gas, the opacity for a given cell can be adjusted based on the evolving dust content. 
This crude approach would allow for a first step investigation using CHYMERA on 
the effects of dust setthng and dust concentration/depletion on the strength of GIs. 

Accretion onto the disk from a surrounding envelope can be simulated, for ex- 
ample, by assuming an Ulrich (1976) flow and by modifying the upper and outer 
boundary cells to permit inflow as well as outflow during fluxing (see Chapter 2). 
However, preliminary work indicates that shocks due to this accretion flow may create 
very hot regions and cause highly disparate radiation and hydrodynamics timescales. 
When radiation algorithms that include some form of cell-to-cell coupling are used 
(M2004, C2006, and the BDNL routines), a numerical catastrophe like the one dis- 
cussed in §2.4.3 could occur. The radiation timescale problem should be addressed 
before envelope accretion is implemented. 

In order to model other accretion mechanisms like the MRI, at least heuristically. 
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an a- viscosity should be included. Eventually, the equations of magnetohydrodynam- 
ics should be added as well. 



8.5.2 Efficiency 

In order to run CHYMERA on Columbia, a NASA Advanced Supercomput- 
ing (NAS) Division computer cluster, a new parallehzation for CHYMERA using 
OpenMP libraries was necessary. Under the old parallehzation, which is the only 
version available at this time through the CVS repository (see CHYMERA's Den), 
the maximum speed-up is about 3-5 x the serial version, depending on system archi- 
tecture, with no gain beyond about eight processors. Because this version of paral- 
lehzation strongly inhibits running CHYMERA in parallel on Columbia, a team of 
specialists from Parallel Software Products, Computer Sciences Corporation, and the 
NAS Division of NASA Ames used a high-end optimization tools on CHYMERA to 
identify inefficiencies in the parallehzation. Using this tool, the team of specialists 
produced a parallel version of CHYMERA that can scale up to 128 processors with a 
speed-up of about 50 x that of the serial version. The new parallehzation is currently 
being added to the repository. 

This underscores the need for periodic efficiency checks of the code. A set of test 
problems and timing tools should be made available to principal CHYMERA users, 
and the results of the efficiency tests should be posted on CHYMERA's Den. It 
should be noted that a message passing version of CHYMERA (distributed mem- 
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ory parallelization) may be released in the near future (see PSP Final Consultancy 
Report on CHYMERA's Den), but there is no guarantee. If a major reorganization 
of CHYMERA occurs within the next several years, the new algorithms should be 
designed with message passing in mind to make a parallelized version of CHYMERA 
portable to distributed as well as shared memory machines. 

There are several ways that the CPU cost of running CHYMERA can be signif- 
icantly reduced; I suggest two. Roughly half of the cells at any time in the code 
contain background densities. At the beginning of every or every few time steps, the 
location of the active grid boundary can be evaluated. If some reliable criterion can 
be determined, then a floating grid boundary could be introduced. Such an algorithm 
could potentially cut the CPU cost of a simulation in half. 

The second suggestion is to implement a multiple time step algorithm. The 
Courant time is almost always determined near the inner disk boundary for these 
simulations, which can be several factors smaller than an appropriate Courant time 
for the middle of the disk. The result is that the mid and outer disk is over-resolved 
temporarily. One can take advantage of these disparate evolution timescales by in- 
troducing a multi-stepping algorithm that advances the inner region of the disk more 
often than the outer regions, similar to what is done for high and low resolution blocks 
in nested grid simulations. 
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8.6 Future Work 

The results of these studies ehcit additional questions, and provide a base for 
future work. 

• What is the effect of complete, 3D cell-to-cell radiative coupling on GIs? 

• How does stellar irradiation affect the effective a when the incident irradiation 
depends in part on the geometry of the disk? 

• What happens to spiral wave dynamics in unstable disks that have both self- 
shadowed regions and regions with normal incident irradiation? 

• How does the effective a in disks vary with star and disk mass? 

• Are the Fourier amplitude spectra and the effective as in these disks converged 
with the resolutions employed? 

• Could dust settling and grain growth or a sudden increase in the EP ionization 

rate due to a nearby source, e.g., a supernova, change the ortho/para statistics? 
What are the dynamic consequences for this change? 

In addition to these questions, the hypothesis that the FU Ori phenomenon and 
chondrule formation are driven by GI activation in dead zones should be explored in 
greater detail. In particular, the slow build up of mass near 1 AU should be followed 
through a burst, and accretion rates as well as shock strengths should be determined. 
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More generally, shock strengths and frequencies should be evaluated for disks in 
asymptotic phases as well as for less violent bursts than explored here. Finally, GIs 
in disks could have cosmochemical consequences, some of which may be observable. 
Gathering fluid element histories for several stages of disk evolution would indicate 
long- and short-term thermodynamic variations that disk material may experience. 



8.7 Concluding Remarks 

Why study protoplanetary disks? One reason is curiosity, a desire to understand 
not just how planetary systems are built, but how the Solar System formed and how 
our civilization developed. A picture of the formation of the Solar System will not 
come from describing any one process. Our story can only be pieced together by 
understanding how multiple astrophysical phenomena are related. Such an investi- 
gation requires the combination of meteoritics, observations, theory, and numerical 
simulations. Throughout history, astronomy has helped shape the way civilizations 
perceive their surroundings, and with the other sciences, has helped us understand 
our place in the universe. It is my sincere hope that this dissertation aids in at least 
some small way to sharpen our cosmic perspective. 



BIBLIOGRAPHY 

Adams, F. C, Lada, C. J., & Shu, F., H. 1987, 312, 788 
Akima, H. 1970, JACM, 17, 4 

Amelin, Y., Krot, A. N., Hutcheon, I. D., & Ulyanov, A. A. 2002, Science, 297, 1678 

Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122 

Andrews, S. M., & WiUiams, J. P. 2005, 619, 175 

Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 

-. 1998, RvMP, 70, 1 

Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 512, 650 

Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213 
Bate, M. R., Ogilvie, G. L, Lubow, S. H., & Pringle, J. E. 2002, MNRAS, 341, 213 
Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924 
Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 

Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton: Princeton 
Univ. Press, 1987) 

Bizzarro, M., Baker, J. A., Haack, & H. 2004, Nature, 431, 275 

263 



264 BIBLIOGRAPHY 
Black, D. C, & Bodenheimer, P. 1975, ApJ, 199, 619 

Bodenheimer, P., Yorke, H. W., Rozyczka, M., & Tohline, J. E. 1990, ApJ, 355, 651 
Boley, A. C, & Durisen, R. H. 2005, LPI, 1286, 8177 
-. 2006, ApJ, 641, 534 

Boley, A. C, Durisen, R. H., Nordlund, A, & Lord, J. 2007c, ApJ, 665, 1254 

Boley, A. C, Durisen, R. H., & Pickett, M. K. 2005, in ASP Conf. Ser. 341, Chondrites 
and the protoplanetary disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San 
Francisco: ASP), 839 

Boley, A. C, Hartquist, T. W., Durisen, R. H., & Michael, S. 2007a, ApJ, 656, L89 

-. 2007b, ApJ, 656, L89 

Boley, A. C, Meji'a, A. C, Durisen, R. H., Cai, K., Pickett, M. K., & D'Alessio, 
P. 2006, 651, 517 

Boss, A. P. 1984a, ApJ, 277, 768 

-. 1984b, MNRAS, 209, 543 

-. 1989, ApJ, 346, 336 

-. 1997, Science, 276, 1836 

-. 1998, ApJ, 503, 923 

-. 2001, ApJ, 562, 367 

-. 2002, ApJ, 567, L149 

-. 2004a, ApJ, 610, 456 



BIBLIOGRAPHY 265 

-. 2004b, ApJ, 616, 1265 
-. 2005, ApJ, 629, 535 
-. 2007, ApJ, 661, L73 

Boss, A. P., & Durisen, R. H. 2005a, ApJ, 621, L137 

-. 2005b, in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed 
A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 821 

Bryden, G., Chen, X., Lin, D. N. C, Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 
514, 344 

Cai, K. 2006, Ph.D. Thesis, Indiana University 

Cai, K., Durisen, R. H., Boley, A. C, Pickett, M. P., Mejia, A. C. 2007, 
arXiv:0706.4046 

Cai, K., Durisen, R. H., Michael, S., Boley, A. C, Mejia, A. C, Pickett, M. K., & 
D'Alessio, P. 2006, ApJ, 636, L149 

Calvet, N., Briceho, C, Hernandez, J., Hoyer, S., Hartmann, L., Sicilia-Aguilar, A., 
Megeath, S. T., & D'Alessio, P. 2005, AJ, 129, 935 

Cameron, A. G. .W. 1978, M&P, 18, 5 

Cassen, P., & Moosman, A. 1981, Icarus, 48, 353 

Chandrasekhar, S. 1960, Radiative transfer (New York: Dover, 1960) 

Cielsa, F. J., & Hood, L. L. 2002, Icarus, 158, 281 

Cohl, H. S., & Tohline, J. E. 1999, 527, 86 



266 BIBLIOGRAPHY 

Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon 
and Breach) 

Cuzzi, J. N. 2004, Icarus, 168, 484 

Cuzzi, J. N., & Alexander, C. M. O'D. 2006, Nature, 441, 483 

Cuzzi, J. N., Ciesla, F. J., Petaev, M. I., Krot, A. N., Scott, E. R. D., & Weiden- 
schilling, S. J. 2005, in ASP Conf. Ser. 341, Chondrites and the protoplanctary 
disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 732 

Cuzzi, J. N., Hogan, R. C, Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496 

D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321 

D'Alessio, P., Calvet, N., Hartmann, L., Franco-Hernandez, R., & Servm, H. 2006, 
ApJ, 638, 314 

D'Alessio, P., Calvet, N., Hartmann, L., Lizano, S., & Canto, J. 1999, ApJ, 527, 893 
Dalgarno, A., Oppenheimer, M., & Black, J. H. 1973, Nature, 245, 100 

Decamph, W. M., Cameron, A. G. W., Bodenheimer, P., & Black, D. C. 1978, ApJ, 

223, 854 

Desch, S. J. 2004, ApJ, 608, 509 

Desch, S. J., & Connolly, H. C, Jr. 2002, M&PSA, 37, 183 

Desch, S. J., & Cuzzi, J. N. 2000, Icarus, 143, 87 

Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485 

DuUemond, C. P., HoUenbach, D., Kamp, I., & D'Alessio, P. 2007, in Protostars and 
Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson: Univ. Arizona Press), 

555 



BIBLIOGRAPHY 267 

Durisen, R. H., Boss, A., Mayer, L., Nelson, A., Quinn, T., & Rice, K. 2007a, in 
Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson: Univ. 

Arizona Press), 607 

Durisen, R. H., Gingold, R. A., Tohline, J. E., & Boss, A. P. 1986, ApJ, 305, 281 

Durisen, R. H., Hartquist, T. W., & Pickett, M. K. 2007b, arXiv: 0709. 1445 

Durisen, R. H., Mejfa, A. C, & Pickett, B. K. 2003, Recent Research Developments 
in Applied Physics, 1, 173 

Durisen, R. H., Yang, S., Cassen, P., & Stahler, S. W. 1989, ApJ, 345, 959 

Eisner, J. A., & Carpenter, J. M. 2006, ApJ, 641, 1162 

Fleming, T. P., & Stone, J. M. 2003, ApJ, 585, 908 

Flower, D. R., & Watt, G. D. 1984, MNRAS, 209, 25 

Flower, D. R., Pineau Des Forets, G., & Walmsley, C. M. 2006, A&A, 449, 621 

Fouchet, T., Lellouch, E., & Feuchtgruber, H. 2004, Icarus, 161, 127 

Fromang, S., Balbus, S. A., Terquem, C., & De Villiers, J. 2004, ApJ, 616, 364 

Fuente, A., Martin-Pintado, J., Rodriquez-Fernandez, N. J., Rodriguez-Franco, A., 
de Vicente, P., & Kunze, D. 1999, ApJ, 518, 45 

Furlan, E., Hartmann, L, Calvet, N., D'Alessio, P., Franco- Hernandez, R., Forrest, 
W. J., Watson, D. M., Uchida, K. I., Sargent B., Green, J. D., Keller, L. D., & 
Herter, T. L. 2006, ApJS, 165, 568 

Gammie, C. F. 1996, ApJ, 457, 355 

-. 2001, ApJ, 553, 174 



268 



BIBLIOGRAPHY 



Gomez, G. C, & Cox, D. P. 2002, ApJ, 580, 235 

Green, J. D., Hartmann, L., Calvet, N., Watson, D. M., Ibrahimov, M., Furlan, E., 
Sargent, B., & Forrest, W. J. 2006, ApJ, 648, 1099 

Greene, T. P., Wilking, B. A., Andre, P., Young, E. T., & Lada, C. J. 1994, ApJ, 
434, 614 

Hachisu, 1. 1986, ApJS, 61, 479 

Barker, D. E., & Desch, S. J. 2002, ApJ, 565, 109 

Bartmann, L. 1998, Accretion Processes in Star Formation (Cambridge: Cambridge 
Univ. Press) 

-. 2005, in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed A. N. Krot, 
E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 1003 

Bartmann, L., D'Alessio, P., Calvet, N., & MuzeroUe, J. 2006, ApJ, 648, 484 

Bartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207 

Bawley, J. F., Smarr, L. L., & Wilson, J. R. 1984, ApJ, 277, 296 

Bayashi, C, Nakazawa, K., Nakagawa, Y. 1985, in Protostars and Planets II, 
ed. D. C. Black, & M. S. Matthews (Tucson: Univ. Arizona Press), 1100 

Beinemann, T., Dobler, W., Nordlund, A, & Brandenburg, A. 2006, A&A, 448, 731 

Berbig, G. B. 1960, ApJS, 4, 337 

-. Berbig, G. B. 1977, ApJ, 217, 693 

Bewins, R. B., Connolly, B. C, Lofgren, G. E. , Jr., & Libourel, G. 2005, in 
ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed A. N. Krot, 
E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 286 



BIBLIOGRAPHY 



269 



Hewins, R., Jones, R., & Scott, E. 1996, Chondrules and the Protoplanetary Disk 
(Cambridge, UK: Cambridge University) 

Home, J. H., & Baliunas, S. L. 1986, ApJ, 302, 757 

Hubeny, I. 1990, ApJ, 351, 632 

Huss, G. R., Alexander, C. M. O'D, Palme, H., Bland, P. A., & Wasson, T. J. 2005, 
in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed A. N. Krot, 
E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 701 

lida. A., Nakamoto, T., Susa, H., & Nakagawa, Y. 2001, Icarus, 153, 430 

Itoh, S., Rubin, A. E., Kojima, H., Wasson, J. T., Yurimoto, H. 2002, LPI, 33.1490 

Jayawardhana, R., Hartmann, L., & Calvet, N. 2001, ApJ, 548, 310 

Johansen, A., & Klahr, H. 2005, 634, 1353 

Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131 

Jones, R. H., Grossman, J. N., & Rubin, A. E. 2005, in ASP Conf. Ser. 341, Chondrites 
and the protoplanetary disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San 
Francisco: ASP), 251 

Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714 

-. 1991, ApJ, 383, 664 

Kitamura, Y., Momose, M., Yokogawa, S., Kawabe, R., Tamura, M., & Ida, S. 2002, 
ApJ, 581, 357 

Korycansky D. G., & Pringle, J. E. 1995, MNRAS, 272, 618 

Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2006, arXiv:astro- 
ph/0611003 



270 BIBLIOGRAPHY 

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, 656, 656 

Kuiper, G. P. 1951, PNAS, 37, 1 

Lada, C. J. 1987, lAUS, 115, 1 

Lada, C. J., & Wilking, B. A. 1984, 287, 610 

Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mechanics (2nd edition: Butterworth- 
Heinemann) 

Larson, R. B. 1984, MNRAS, 206, 197 

Laughlin, G., & Bodenheimer, P. 1994, 436, 335 

Laughlin, G., & Korchagin, V. 1996, 480, 855 

Laughlin, G., Korchagin, V., & Adams, F. C. 1998, 504, 945 

Laughhn, G., & Rozyczka, M. 1996, ApJ, 456, 279 

Lawson, W. A., Feigelson, E. D., & Huenemoerder, D. P. 1996, MNRAS, 280, 1071 

Le Bourlot, J. 2000, A&A, 360, 656 

Lin, D. N. C., & Papaloizou, J. 1980, MNRAS, 191, 37 

Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990, ApJ, 364, 326 

Lissauer, J. J. 1987, Icarus, 69, 249 

Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630 

Lubow, S. H. 1981, ApJ, 245, 274 

Lubow, S. H., & Ogilvie, G. I. 1998, ApJ, 504, 983 

Lubow, S. H., & Pringle, J. E. 1993, ApJ, 409, 360 



BIBLIOGRAPHY 271 
Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1 

MacPherson, G. J., Simon, S. B., Davis, A. M., Grossman, L., & Krot, A. N. 2005, 
in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed A. N. Krot, 
E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 225 

Martos, M. A., & Cox, D. P. 1998, ApJ, 509, 703 

Massey, B. S. 1970, Mechanics of Fluids (2nd ed.; London: Van Norstrand Reinhold) 
Mayer, L., Lufkin, G., Quinn, T., & Wadsley J. 2007, ApJ, 661, 77 
Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2004, ApJ, 609, 1045 
McKeegan, K. D. et al. 2006, Science, 314, 1724 
Mejfa, A. C. 2004, Ph.D. Thesis, Indiana University 

Mejia, A. C, Durisen, R. H., Pickett, M. K., & Cai, K. 2005, ApJ, 619, 1098 

Mihalas, D., & Weibel-Mihalas, B. 1984, Foundations of radiation hydrodynamics 
(New York: Oxford University Press, 1984) 

Miura, H., & Nakamoto, T. 2006, ApJ, 651, 1272 

Miyake, K., & Nakagawa, Y. 1995, ApJ, 441, 361 

Monaghan, J. J. 1992, ARA&A, 30, 543 

Morfill, G. E., Durisen, R. H., & Turner, G. W. 1998, Icarus, 134, 180 

MuzeroUe, J., Calvet, N., Hartmann, L., & D'Alessio, P. 2003, ApJ, 597, L149 

Natta, A., Testi, L., Calvet, N., Henning, T., Waters, R., & Wilner, D. 2007, in 
Protostars and Planets V, ed. B. Reipurth, D. Jewitt, and K. Keil (Tucson: Univ. 
Arizona Press), 767 



272 BIBLIOGRAPHY 

Nelson, A. F. 2000, ApJ, 537, 65 
-. 2006, MNRAS, 373, 1039 

Nelson, A. F., Benz, W., & Ruzmaikina, T. V. 2000, ApJ, 529, 357 

Norman, M. L., & Winkler, K. H. A. 1986, in NATO Advanced Research Workshop 
on Astrophysical Radiation Hydrodynamics 

Ogilvie, G. I. 1998, MNRAS, 297, 291 

-. 2002a, MNRAS, 330, 937 

-. 2002b, MNRAS, 331, 1053 

Oishi, J. S., Mac Low, M., Menou, K. 2007, arXiv:astro-ph/0702549 

Osorio, M., D'Alessio, R, MuzeroUe, J., Calvet, N., & Hartmann, L. 2003, ApJ, 586, 
1148 

Osterbrock, D. E. 1962, ApJ, 136, 359 

Padgett, D. L., Brandner, W., Stapelfeldt, K. R., Strom, S. E., Tereby, S., & Koerner, 
D. 1999, AJ, 117, 1490 

Pathria, R. K. 1996, Statistical Mechanics (2nd ed.; Oxford: Butterworh-Heinemann) 

Petaev, M. I., & Wood, J. A. 2005, in ASP Conf. Ser. 341, Chondrites and the 
protoplanetary disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: 
ASP), 373 

Palla, F., & Stabler, S. W. 1990, ApJ, 360, 47 
Pickett, B. K. 1995, Ph.D. Thesis, Indiana University 

Pickett, B. K., Cassen, P. M., , Durisen, R. H., & Link, R. P. 1998, ApJ, 504, 468 



BIBLIOGRAPHY 273 
-. 2000a, ApJ, 529, 1034 

Pickett, M. K., & Durisen, R. H. 2007, ApJ, 654, L155 

Pickett, B. K., Durisen, R. H., Cassen, P. M., & Mejia A. C. 2000b, ApJ, 540, L95 

Pickett, B. K., Durisen, R. H., & Davis, G. A. 1996, ApJ, 458, 714 

Pickett, B. K., Durisen, R. H., & Link, R. 1997, Icarus, 126, 243 

Pickett, B. K., Meji'a, A. C, Durisen, R. H., Cassen, P. M., Berry, D. K., & Link, 
R. P. 2003, ApJ, 590, 1060 

Pilipp, W., Hartquist, T. W., MorfiU, G. E., & Levy, E. H. 1998, A&A, 331, 121 

Pollack, J. B., Hollenbach, D., Beckwith, S., Simonelli, D. P., Roush, T., & Fong, 
W. 1994, 421, 615 

Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Numerical recipes. The art 
of scientific computing (Cambridge: University Press) 

Pringle, J. E. 1981, ARA&A, 19, 137 

Rafikov, R. R. 2005, ApJ, 621, L69 

-. 2007, ApJ, 662, 642 

Reipurth, B. 2005, in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, 
ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 54 

Rice, W. K. M, Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 
1025 

Rice, W. K. M, Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56 
Roberts, W. W., Huntley, J. M., & van Albada, G. D. 1979, ApJ, 233, 67 



274 



BIBLIOGRAPHY 



Rodriguez-Fernandez, N. J., Martm-Pintado, J., de Vicente, P., Puente, A., 
Hiittemeister, S., Wilson, T. L., & Kunze, D. 2000, A&A, 356, 695 

Ruden, S. P., & Pollack, J. B. 1991, ApJ, 375, 740 

Russell, S. S., Krot, A. N., Huss, G. R., Keil, K., Itoh, S., Yurimoto, H., & Macpher- 
son, G. J. 2005, in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, 
ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 317 

Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 

Scott, E. R. D., & Krot, A. N. 2005, in ASP Conf. Ser. 341, Chondrites and the 
protoplanetary disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco: 
ASP) 

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 

Shu, F. H., Shang, H., Gounelle, M., Glassgold, A. E., & Lee, T. 2001, ApJ, 548, 
1029 

Sod, G. A. 1978, JCoPh, 27, 1 

Spitzer, L., & Tomasko, M. G. 1968, ApJ, 152, 971 

Stahler, S. W. 1983, ApJ, 274, 822 

Stahler, S. W. 1988, ApJ, 332, 804 

Stamatellos, D., Whitworth, A. P., & Ward-Thompson, D. 2007, MNRAS, 379, 1390 

Stepinski, T. F. 1992, Icarus, 97, 130 

Sternberg, A., & Neufeld, D. A. 1999, ApJ, 516, 371 

Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 



BIBLIOGRAPHY 275 

Strom, S. E., & Edwards, S. 1993, ASPC, 36, 235 

Tachibana, S., & Huss, G. R. 2005, GeCoA, 69, 3075 

The, P. S., Perez, M. R., & van den Heuvel, E. P. J. 1994, ASPC, 62, 23 

Tohline, J. E. 1980, ApJ, 235, 866 

Tomley, L., Cassen, P., & Steiman-Cameron, T. 1991, ApJ, 382, 530 

Tomley, L., Steiman-Cameron, T. Y., Cassen, P. 1994, ApJ, 422, 850 

Toomre, A. 1964, ApJ, 139, 1217 

Ulrich, R. K. 1976, ApJ, 210, 377 

Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617 

van Albada, G. D., van Leer, B., & Roberts, W. W. 1982, A&A, 108, 76 

van den Ancker, M. E., The, P. S., Tjin A Djie, H. R. E., Catala, C, de Winter, D., 
Blondel, P. F. C, Waters, L. B. F. M. 1997, A&A, 324, 33 

Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, 137 

-. 2006, ApJ, 650, 956 

Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astronomy, 9, 137 

Walmsley, C. M., Flower, D. R., & Pineau des Forets, G. 2004, A&A, 418, 1035 

Webber, W. R. 1998, ApJ, 506, 329 

WeidenschiUing, S. J. 1977, Ap&SS, 51, 153 

WiUiams, H. A. 1988, Ph.D. Thesis, Louisiana State University 

Whipple, F. L. 1966, Science, 153, 54 



276 



BIBLIOGRAPHY 



Whitehouse, S. C, & Bate, M. R. 2006, MNRAS, 367, 32 

Wood, J. A. 1963, Icarus, 2, 152 

-. 1996, Meteoritics Planet. Sci., 31, 641 

-. 2005, in ASP Conf. Ser. 341, Chondrites and the protoplanetary disk, ed A. N. Krot, 
E. R. D. Scott, & B. Reipurth (San Francisco: ASP), 953 

Wooden, D., Barker, D. E., & Brearley, A. J. 2005, in ASP Conf. Ser. 341, Chondrites 
and the protoplanetary disk, ed A. N. Krot, E. R. D. Scott, & B. Reipurth (San 
Francisco: ASP), 774 

Yang, S. X. 1992, Ph.D. Thesis, Indiana University 

Yorke, H. W., Bodenheimer, P., & Laughhn, G. 1993, ApJ, 411, 274 

Zhu, Z., Hartmann, L., Calvet, N., Hernandez, J., MuzeroUe, J., & Tannirkulam, 
A. 2007, arXiv:0707.3429 



Aaron Christopher Boley 



CONTACT INFORMATION 



Swain Hall West 319 
Department of Astronomy 
Indiana University 
Bloomington, IN 
47405-7105, USA 



Phone: 812.855.6911 

Fax: 812.855.8725 

E-mail: acboley ©astro . indiana.edu 

Web site: http://orion.astro.indiana.edu/~aaron 



Education 



Indiana University, 

Bloomington, IN, USA 

Ph.D. Astrophysics, 
2007 

- Advisor: Richard H. Durisen 



Mount Union College, 

AUiance, OH, USA 

September B.S. Physics & Astronomy, German, May 

2002 

- Magna Cum Laude 
M.A. Astronomy, October 2004 ' Minor Mathematics 



Research Interests 

Computational/theoretical astrophysics: chondrule formation and dust processing, planet 
formation, protoplanetary disk evolution, and radiative hydrodynamics. 

Research Experience 



Indiana University: (Supervisor - Richard H. Durisen) 

Graduate student, Ph.D. research 
NASA GRC, Cleveland, OH: (Supervisor - Jeffrey Wilson) 

Computational Modeling of Left-Handed Metamaterials 
NRAO, Socorro, NM: (Supervisor - Mark Claussen) 

Methanol Masers in Star Forming Regions 

Honors and Awards 



2002-2007 



Summer 2002 



Summer 2001 



• NASA Graduate Student Research Program Fellow, 2005-2007 

• Indiana University Astronomy Department's Research Recognition, 2006 

• Indiana Space Grant Fellow, 2005 

• Indiana University Astronomy Department's Teaching Assistant Recognition, 2004 



Teaching Experience 



• NSF Sponsored Research Experience for Undergraduates (REU) Mentor, June-July 
2005 & 2006 

• Associate Instructor: AlOO Introduction to the Solar System and A105 Stars and 
Galaxies, 6 Sessions between 2003 and 2005 

• Teaching Assistant: AlOO and General Astronomy II, Sept.-Dec. 2002, Jan.-May 2003 

Public Outreach Activities 

• Kirkwood Open House, host for public observing sessions, 2002-present 

• Physics and Astronomy Open House, hosted astronomy activities for the public, Fall 
2002-2006 

• Brownie Math and Science Day, discuss various astronomy topics with children. Fall 
2003 and 2004 

• Science Olympiad, test writer/grader, Spring 2002-2005 

Refereed Publications 

8 "Gravitational Instabilities, Chondrule Formation, and the FU Orionis Phenomenon," 
A. C. Boley & R. H. Durisen 2007, ApJ Letters, submitted. 

7 "The Thermal Regulation of Gravitational Instabilities in Protoplanetary Disks. IV. 

Simulations with Envelope Irradiation," 

K. Cai, R. H. Durisen, A. C. Boley, M. K. Pickett, k A. C. Mejfa 2007, ApJ, 
submitted. 

6 "3D Radiative Hydrodynamics for Disk Stability Simulations: A Proposed Testing 
Standard," 

A. C. Boley, R. H. Durisen, A. Nordlund, & J. Lord 2007, ApJ, 665, 1254. 

5 "The Internal Energy for Molecular Hydrogen in Gravitationally Unstable Protoplan- 
etary Disks," 

A. C. Boley, T. W. Hartquist, R. H. Durisen, & S. Michael 2006, ApJ, 656, L89. 

4 "The Thermal Regulation of Gravitational Instabilities in Protoplanetary Disks III. 

Simulations with Radiative Cooling & Realistic Opacities," 

A. C. Boley, A. C. Mejfa, R. H. Durisen, M. K. Pickett, k P. D'Alessio 2006), ApJ, 
651, 517. 

3 "Hydraulic/Shock- Jumps in Protoplanetary Disks," 
A. C. Boley k R. H. Durisen 2006, ApJ, 641, 534. 



2 "The Effects of Metallicity &; Grain Size on Gravitational Instabilities in Protoplan- 

etary Disks," 

K. Cai, R. H. Durisen, S. Michael, A. C. Boley, A. C. Mejia, M. K. Pickett, & 
P. D'Alessio 2006, ApJL, 636, 149. 

1 "The Three-Dimensionality of Spiral Shocks in Disks: Did Chondrules Catch a Break- 
ing Wave?" 

A. C. Boley, R. H. Durisen, k. M. K. Pickett 2005, In Chondrites and the Proto- 
planetary Disk, ASPC Series, 341, 839. 

Selected Conference Presentations 

7 "A Test Suite for 3D Radiative Hydrodynamics Simulations of Protoplanetary Disks," 
A. C. Boley, R. H. Durisen, A. Nordlund, & J. Lord 2006, BAAS, 38, 995. 

6 "3D Radiative Hydrodynamics Simulations of Protoplanetary Disks: A Comparison 
Between Two Radiative Cooling Algorithms," 
J. Lord, A. C. Boley, & R. H. Durisen 2006, BAAS, 38, 996 

5 "The Effects of Varied Initial Conditions on Protoplanetary Disks," 
S. Michael, A. C. Boley, & R. H. Durisen 2006, BAAS, 38, 1052. 

4 "Hydraulic/Shock- Jumps in Protoplanetary Disks," 
A. C. Boley & R. H. Durisen 2005, BAAS, 37, 1164. 

3 "Star Formation in the Outer Disk of NGC 5964," 

C. P. McKinney, A. C. Boley, L. van Zee, D. Schade, & S. Cote 2005, BAAS, 37, 
1392. 

2 "Linking Chondrules to the Formation of Jupiter Through Nebular Shocks," 
A. C. Boley & R. H. Durisen 2005, LPI 1286, 8177. 

1 "The Three-Dimensionality of Spiral Shocks in Disks: Did Chondrules Catch a Break- 
ing Wave?" 

A. C. Boley, R. H. Durisen, M. K, Pickett 2004, Workshop on Chondrites & the 
Protoplanetary Disk, 9016. 



